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Preface 


The  computer  program  described  in  this  report  evolved  during  the  course  of  a 
number  of  projects  in  the  Semiconductor  Technology  Program  in  the  Electron 
Devices  Division  at  the  National  Bureau  of  Standards.     This  program  serves  to 
focus  NBS  efforts  to  enhance  the  performance,   interchangeability,  and 
reliability  of  discrete  semiconductor  devices  and  integrated  circuits  through 
improvements  in  measurement  technology  for  use  in  specifying  materials  and 
devices  in  national  and  international  commerce  for  use  by  industry  in 
controlling  device  fabrication  processes.     Its  major  thrusts  are  the 
development  of  carefully  evaluated  and  well-doctmiented  test  procedures  and 
associated  technology  and  the  dissemination  of  such  information  to  the 
electronics  community. 
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Semiaonductor  Measurement  Teahnology: 
A  FORTRAN  Program  for  Calculating  the  Electrical  Parameters 

of  Extrinsic  Silicon 

by 

R.  D.  Larrabee,  W.  R.  Thurber,   and  W.  M.  Bullis 
National  Bureau  of  Standards 
Washington,   DC  20234 

ABSTRACT 

Many  electrical  properties  of  silicon  are  strongly  depen- 
dent upon  the  specific  nature  and  density  of  the  active 
impurities  present.     Calculation  of  these  electrical 
properties  hinges  on  the  solution  of  the  charge  balance 
equation  to  determine  the  position  of  the  Fermi  level  for 
each  specific  case  of  interest.     A  FORTRAN  program  is 
presented  that  performs  this  determination  and  then  cal- 
culates some  of  the  often-used  electrical  parameters  of 
silicon  as  a  function  of  temperature.     Results  obtained 
from  this  program  have  proven  useful  in  interpreting  Hall 
effect  data,   determining  the  degree  of  ionization  of  the 
separate  dopant  states  as  a  function  of  temperature,  pre- 
dicting the  behavior  of  specimens  when  the  dopant  picture 
is  intentionally  (or  conceptually)  changed,   and  under- 
standing the  variations  in  the  relative  roles  of  the  dif- 
ferent scattering  mechanisms  on  carrier  mobility  as  the 
temperature  is  changed. 

Key  Words:     Carrier  density;  computer  program;  electri-  - 
cal  properties  of  silicon;   Hall  effect;  mobility;  resis- 
tivity; silicon. 

1 .  INTRODUCTION 

Trace  impurities  can  act  as  donors  or  acceptors  and  thus  contribute  mobile 
charge  carriers  that  can  have  a  profound  effect  on  resistivity  and  other 
electrical  properties.     Since  carrier  density  and  mobility  are  functions  of 
temperature  as  well  as  the  density  and  type  of  impurities  present,   it  is  im- 
practical to  provide  general  tables  listing  all  possible  situations  that 
might  occur  in  practice.     Therefore,   it  becomes  necessary  to  calculate  the 
desired  electrical  properties  for  each  impurity  situation  and  temperature  of 
interest. 

The  FORTRAN  program  presented  in  this  report  was  written  to  accomplish  this 
objective  for  the  case  of  simple  donor  and  acceptor  states  in  silicon  (here- 
after called  dopant  states).     The  density,   activation  energy,   and  degeneracy 
of  all  active  donors  and  acceptors  are  initially  provided  as  inputs  to  the 
program.     These  data  are  then  used  by  the  program  to  calculate  the  position 
of  the  Fermi  level  by  solving  the  detailed  charge  balance  equation  (i.e.,  the 
sum  of  the  positive  charges  in  the  valence  band  and  ionized  donors  equals  the 
sum  of  negative  charges  in  the  conduction  band  and  ionized  acceptors).  After 


the  position  of  the  Fermi  level  has  been  found,   the  program  computes  many  of 
the  more  often  used  electrical  properties  of  silicon  that  depend  on  carrier 
density  or  carrier  scattering  (e.g.,   resistivity,   carrier  mobility,   and  Hall 
coefficient).     The  scattering  mechanisms  which  are  included  in  the  calcula- 
tion of  the  carrier  mobility  are  lattice,   ionized  impurity,   neutral  impurity, 
and  electron-electron  (or  hole-hole)  interactions.     Both  electron  and  hole 
mobilities  are  calculated  and  an  effective  conductivity  mobility  is  calcu- 
lated as  the  sum  of  the  electron  and  hole  mobilities  weighted  by  the  density 
of  each.     The  effective  conductivity  mobility  computed  this  way  is  the  same 
as  the  majority  carrier  conductivity  mobility  in  extrinsic  specimens.     In  in- 
trinsic specimens,   it  corresponds  to  a  hypothetical  mobility  for  the  total 
carrier  density  (i.e.,   electron  plus  hole  density).     The  expressions  for  cal- 
culating the  resistivity  and  Hall  coefficient  include  appropriately  weighted 
contributions  from  both  charge  carriers.     In  the  absence  of  agreement  as  to 
appropriate  values,   the  Hall  scattering  factors  (i.e.,   ratio  of  Hall  to  con- 
ductivity mobility)  for  both  electrons  and  holes  were  taken  to  be  unity. 

In  its  present  form,   the  program  handles  up  to  five  independent  impurity 
states,  two  donors  and  three  acceptors.     The  donor  levels  are  assumed  to  be 
neutral  when  occupied  by  an  electron  and  to  have  a  temperature-independent 
activation  energy  referenced  to  the  conduction  band  edge.     The  acceptor  lev- 
els are  assumed  to  be  neutral  when  occupied  by  a  hole  and  to  have  a 
temperature-independent  activation  energy  referenced  to  the  valence  band 
edge.     Provision  is  made  for  the  reduction  in  activation  energy  for  both 
donors  and  acceptors  at  the  higher  dopant  densities.     These  assumptions  are 
appropriate  for  impurities  from  column  3  and  5  of  the  periodic  table  (e.g., 
boron  and  phosphorus)  and  may,   or  may  not,   be  appropriate  for  other  impuri- 
ties with  larger  activation  energies  or  multiple  energy  levels. 

The  key  algorithms  in  this  program  are  the  routine  for  calculating  the  posi- 
tion of  the  Fermi  level  and  the  routine  for  calculating  the  carrier  mobili- 
ties.    These  algorithms  are  most  accurate  over  the  temperature  interval  from 
100  to  500  K  and  for  total  doping  densities  below  10^^  cm~^.  These 
two  algorithms  have  been  written  as  subroutines  in  the  present  program  so 
that  they  can  be  used  as  a  starting  point  in  the  development  of  other  pro- 
grams that  require  these  algorithms.     Notice,   however,   that  the  mobility  rou- 
tine requires  knowledge  of  many  of  the  parameters  calculated  by  the  Fermi 
level  routine,  and  thus  it  cannot  be  used  without  first  executing  the  Fermi 
level  routine. 

Variations  of  the  present  program  have  proven  useful  for  interpreting  Hall 
data   [1],  predicting  the  electrical  behavior  when  the  dopant  picture  (e.g., 
compensation)  is  changed  [2,3],   determining  the  degree  of  ionization  of  the 
active  dopant  state  in  extrinsic  infrared  detectors  [4] ,   and  understanding 
the  changing  role  of  the  various  scattering  mechanisms  on  mobility  as  the 
temperature  or  dopant  density  changes. 

The  following  section  provides  the  information  necessary  to  use  the  program. 
Section  3  describes  the  various  equations  used  in  the  program  so  that  the  in- 
terested user  may  verify  that  they  are  appropriate  for  his  particular  appli- 
cation.    A  detailed  listing  of  the  complete  program  is  given  in  section  4. 
This  listing  is  heavily  annotated  so  that  the  user  can  identify  each  equation 
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and  follow  each  step  of  the  program.     Section  5  provides  examples  of  typical 
output  formats. 

2.      DESCRIPTION  OF  PROGRAM  OPERATION 
The  present  program  can  be  divided  into  four  main  parts  as  follows: 

1.  Read  data  cards  and  set  up  options. 

2.  Compute  position  of  the  Fermi  level,   occupancy  of  all  dopant  states, 
and  density  of  carriers  in  the  conduction  and  valence  bands. 

3.  Compute  carrier  mobilities  and  related  parameters. 

4.  Print  output  listings  and  reinitialize  for  next  run. 

In  its  present  form,   the  program  can  handle  up  to  two  donor  states  (called  Dl 
and  D2)  and  up  to  three  acceptor  states  (called  Al,   A2,   and  A3).     For  each 
dopant  state,   the  program  expects  to  find  a  data  card  specifying: 

1.  the  name  of  the  dopant  state  (up  to  12  alphanumeric  characters), 

2.  the  density  of  the  dopant  state  in  reciprocal  cubic  centimeters, 

3.  the  activation  energy  of  the  dopant  state  at  infinite  dilution  in 
electron  volts,*  and 

4.  the  degeneracy  of  the  dopant  state  (an  integex-  >  unity). 

The  format  of  these  data  cards  is  2A6,   E10.4,   2F10.4,   corresponding  to  the 
four  items  of  input  data  -  name,  density,   activation  energy,   and  degeneracy, 
respectively.     The  format  for  the  name   (2A6)  assumes  that  the  machine  can 
store  12  alphanumeric  characters  in  two  6-character  words.     This  may  have  to 
be  changed  for  machines  with  different  alphanumeric  storage  formats.     If  a 
particular  dopant  state  is  absent,   its  data  card  must  still  be  present  to 
inform  the  program  of  this  fact.     The  data  cards  representing  missing  dopant 
states  use  the  name  DUMMY  and  the  rest  of  the  card  is  left  blank.     The  order 
of  the  five  data  cards  is: 


1. 

Dl 

-  first  donor  state. 

2. 

D2 

-  second  donor  state. 

3. 

Al 

-  first  acceptor  state. 

4. 

A2 

-  second  acceptor  state,  and 

5. 

A3 

-  third  acceptor  state. 

It  is  recommended  that  this  same  order  (i.e.,  Dl  through  A3)  be  the  order  of 
occurrence  of  the  corresponding  dopant  states  from  the  conduction  band  to  the 
valence  band  (if  possible)  to  aid  in  remembering  which  symbol  (i.e.,  Dl 
through  A3)  represents  each  dopant  state. 


*For  donor  states,   the  activation  energy  is  the  magnitude  of  the  energy 
difference  between  the  conduction  band  and  the  dopant  state;  for  acceptor 
states,   it  is  the  magnitude  of  the  energy  difference  between  the  valence 
band  edge  and  the  dopant  state.     Both  of  these  quantities  are  assijmed  to  be 
independent  of  temperature. 
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In  the  normal  (no-option)  output  format,  the  program  lists  the  following 
quantities  on  the  output  page  for  each  set  of  five  input  data  cards  sup- 
plied: 

1 .  the  position  of  the  Fermi  level  with  respect  to  the  valence  band  in 
electron  volts, 

2.  the  density  of  electrons  in  the  conduction  band  in  reciprocal  cubic 
centimeters, 

3.  the  density  of  holes  in  the  valence  band  in  reciprocal  cubic 
centimeters, 

4.  the  effective  conductivity  mobility  in  square  centimeters  per 
volt-second, 

5.  the  Hall  coefficient  in  cubic  centimeters  per  coulomb, 

6.  the  electrical  resistivity  in  ohm-centimeters,  and 

7.  the  ionized  density  of  dopant  species  Dl  through  A2  in  reciprocal 
cubic  centimeters.     (The  DOUBLE  ON  option  discussed  below  can  be 
used  to  output  the  ionized  density  of  species  A3. ) 

These  ten  items  are  printed  across  a  single  line  on  the  output  listing  in  the 
order  given.     Since  these  parameters  are  functions  of  temperature,   the  output 
listing  consists  of  a  series  of  such  lines,   each  preceded  by  its  temperature 
and  the  corresponding  value  of  1000/T.     The  resulting  table  is  provided  with 
descriptive  column  headings  and  documentation  of  the  appropriate  input  infor- 
mation from  the  data  cards  (see  first  sample  output  listing  in  sec.  5).  The 
program  is  capable  of  listing  additional  output  information,   as  discussed  in 
the  section  on  options  to  follow. 

The  above  quantities  are  evaluated  and  listed  as  a  function  of  temperature 
over  the  range  from  4.2  to  500  K.     This  temperature  range  is  divided  into 
three  regions,  with  break-points  at  77  K  (liquid  nitrogen  temperature)  and 
300  K  (room  temperature).     The  program  will  perform  its  evaluations  at  each 
break-point  (i.e.,   77  and  300  K) ,   as  well  as  at  the  extreme  points  (i.e.,  4.2 
and  500  K) .     After  completion  of  the  first  evaluation  at  4.2  K,   the  program 
jumps  to  some  higher  specified  temperature  and  then  starts  to  increase  the 
temperature  with  equal  increments  of  1 000/( temperature ) .     The  incremental 
values  of  1000 /(temperature)  are  specified  independently  in  each  of  the  three 
temperature  regions.     As  presently  written,   the  temperature  following  4.2  K 
is  20  K  (1000/T  =  50),  and  then  the  temperature  is  increased  in  intervals  of 
1000/T  =  2  until  the  first  break-point  at  77  K  is  reached.     Within  the  inter- 
val from  77  K  to  300  K,   the  temperature  is  increased  in  intervals  of  1000 A  = 
1.     Above  300  K,   the  temperature  is  increased  in  intervals  of  1000/T  =  0.5. 
The  FORTRAN  parameter  RUNNO  keeps  track  of  things  as  follows: 

RUNNO  =  1,  first  time  through,   T  =  4.2  K, 

RUNNO  =  2,  first  region  up  to  77  K, 

RUNNO  =  3,  second  region  from  77  to  300  K,  and 

RUNNO  =  4,  third  region  above  300  K. 

Many  of  these  details  of  the  temperature  sequence  can  easily  be  changed  by 
simply  changing  the  appropriate  defining  statements  in  the  program  (e.g.,  the 
first  temperature  of  4.2  K  is  defined  by  the  FORTRAN  statement  numbered  170 
in  the  program  listing  in  sec.  4.). 
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2.1  OPTIONS 


Several  options  are  built  into  this  program  to  allow  one  to  modify  the  output 
listing  without  changing  the  program  pev  se.     These  options  are  exercised 
when  a  control  card  is  inserted  into  the  data  card  deck  so  that  it  is  read  by 
the  program  as  a  D1  data  card.     The  program  recognizes  that  it  is  a  control 
card  and  not  a  D1  data  card  by  the  name.     When  such  a  control  card  is  encoun- 
tered,  the  program  sets  up  the  option  desired  and  then  goes  back  and  expects 
the  next  card  to  be  a  D1  data  card.     If  no  control  cards  are  included  in  the 
deck  of  input  data  cards,   the  program  simply  generates  an  output  listing  as 
outlined  previously. 

This  is  considered  to  be  the  normal  (or  no-option)  output  format  and  is  the 
format  that  should  suffice  for  most  applications  of  this  program. 

The  control  cards  consist  of  a  single  (or  double)  word  name  selected  from  a 
reserved  name  list  designed  not  to  conflict  with  dopant  state  names.     Some  of 
the  control  cards  additionally  contain  numerical  arguments  that  are  needed  to 
specify  the  parameters  of  the  option  desired.     In  either  event,   the  format  of 
all  control  cards  is  identical  to  that  of  the  dopant  state  data  cards  (i.e., 
2A6,  ElO.4,   2F10.4),   so  the  names  and  any  numerical  arguments  must  be  placed 
on  the  control  card  accordingly. 

2.1.1     DOUBLE  ON 

This  option  doubles  the  amount  of  data  printed  on  the  output  listing  at 
each  temperature.     The  format  of  the  output  listing  consists  of  two 
lines  of  output  information  at  each  temperature  instead  of  just  one 
(thus  the  name  DOUBLE  ON).     The  first  line  of  information  is  identi- 
cal to  a  no-option  run*  and  the  second  line  contains  the  parameters 
listed  below.     The  corresponding  control  card  consists  of  the  two  words 
DOUBLE  ON  in  the  first  nine  positions  with  the  rest  of  the  card  blank. 
This  control  card  turns  on  this  option  for  all  subsequent  sets  of  dopant 
state  data  cards,   and  it  remains  on  until  a  DOUBLE  OFF  control  card  is 
read.     This  means  that  the  DOUBLE  ON  control  card  must  be  the  first  card 
of  a  group  of  data  card  sets  for  which  the  double  output  is  desired. 
The  DOUBLE  OFF  control  card  must  be  the  card  which  follows  the  last  data 
card  in  the  group  (i.e.,   the  first  card  of  a  run  for  which  single  output 
only  is  desired) .     In  addition  to  all  the  parameters  listed  in  a  no- 
option  run,   the  DOUBLE  ON  option  will  cause  the  program  to  list  the  fol- 
lowing parameters: 

1.  The  energy  gap  of  silicon  in  electron  volts. 

2.  A  parameter  called  TEST  which  indicates  how  close  the  calculated 
Fermi  level  position  comes  to  solving  the  detailed  charge  balance 
equation.     TEST  is  the  number  of  residual  charges  per  cubic 
centimeter, 

3.  The  voltage  equivalent  of  the  temperature  in  volts  (i.e.,   kT/q) . 

4.  The  position  of  the  Fermi  level  with  respect  to  the  conduction  band 
in  electron  volts. 


*A  run  is  herein  defined  as  the  evaluation  of  one  set  of  five  dopant  state 
data  cards  over  the  full  temperature  range. 
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5.  The  ionized  density  of  the  dopant  state  A3  in  reciprocal  cubic 
centimeters  (there  was  insufficient  room  for  this  to  be  included  on 
the  single-line  no-option  listing) . 

6.  The  electrical  conductivity  in  mho  per  centimeter. 

7.  The  mobility  of  electrons  in  the  conduction  band  in  square  centime- 
ters per  volt-second,   as  follows: 

A.  the  upper  limiting  mobility  of  electrons  set  by  lattice  and 
electron-electron  scattering, 

B.  the  mobility  obtained  by  adding  the  effects  of  ionized  impurity 
scattering  (corrected  for  electron-electron  scattering)  to  A, 
above,  and 

C.  the  mobility  obtained  by  adding  the  effects  of  neutral  impurity 
scattering  to  B,  above. 

8.  The  mobility  of  holes  in  the  valence  band  in  square  centimeters  per 
volt-second  as  follows: 

A.  the  upper  limiting  mobility  of  holes  set  by  lattice  and  hole- 
hole  scattering, 

B.  the  mobility  obtained  by  adding  the  effects  of  ionized  impurity 
scattering  (corrected  for  hole-hole  scattering)  to  A,  above, 
and 

C.  the  mobility  obtained  by  adding  the  effects  of  neutral  impurity 
scattering  to  B,  above. 

These  parameters  were  selected  to  complement  the  parameters  of  a  no- 
option  run  and  are  listed  on  the  line  immediately  below  them  on  the  out- 
put listing,   as  shown  in  section  5.3.     The  sets  of  two-line  data  corre- 
sponding to  each  temperature  are  separated  by  two  blank  lines  in  order 
to  make  it  more  obvious  which  of  the  two  coliamn  headings  applies  to  each 
line.     Some  of  the  parameters  listed  here  are  also  listed  in  the  CON- 
STANTS option  discussed  below,   and  some  are  modified  versions  of  parame- 
ters listed  in  the  no-option  listing  (e.g.,   the  position  of  the  Feirmi 
level  with  respect  to  the  conduction  band  instead  of  with  respect  to  the 
valence  band).     The  three  mobility  values  listed  for  electrons  and  holes 
provide  information  needed  for  giving  one  a  picture  of  the  changing  role 
of  the  various  scattering  mechanisms  with  temperature.     Although  the  ma- 
jority carrier  mobility  is  calculated  as  accurately  as  possible,  less 
attention  has  been  paid  to  the  calculation  of  minority  carrier  mobility. 
(See  sec.  3.2  below.) 

2.1.2  DOUBLE  OFF 

This  control  card  turns  off  the  DOUBLE  ON  option  as  explained  above. 
The  control  cards  DOUBLE  ON  and  DOUBLE  OFF  simply  set  and  clear  a 
switch,   respectively,   and  affect  nothing  but  the  format  of  the  output 
listing.     Therefore,   it  is  permissible  to  use  any  of  the  other  control 
card  options  while  in  the  DOUBLE  ON  mode  of  output  listing.     The  default 
condition  (i.e.,   if  neither  DOUBLE  ON  nor  DOUBLE  OFF  control  cards 
appear  in  the  deck  of  data  cards)   is  DOUBLE  OFF. 

2.1.3  CONSTANTS 

This  option  generates  a  listing  of  some  of  the  parameters  used  in  the 
evaluation  that  are  not  functions  of  dopant  density.     These  parameters 
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are  evaluated  at  the  same  temperatures  as  a  no-option  run.     The  corre- 
sponding control  card  consists  of  the  single  word  CONSTANTS  in  the  first 
nine  positions  with  the  rest  of  the  card  blank.     The  CONSTANTS  option 
does  not  relate  in  any  way  to  any  dopant  state  data  cards  that  might 
follow.     Consequently,   the  next  card  could  be  another  control  card  or 
the  D1  data  card  of  another  run.     The  CONSTANTS  control  card  can  be  lo- 
cated anywhere  in  the  deck  of  cards  to  be  read  by  the  program  as  long  as 
it  is  encountered  when  the  program  expects  to  find  a  D1  data  card.  The 
CONSTANTS  control  card  causes  the  program  to  list  the  following 
parameters: 

1.  100 0/( temperature )  in  inverse  kelvin, 

2.  the  temperature  in  kelvin, 

3.  the  voltage  equivalent  of  the  temperature  in  volts  (i.e.,  kT/q), 

4.  the  band  gap  in  silicon  in  electron  volts, 

5.  the  relative  dielectric  constant  of  silicon, 

6.  the  upper  limiting  mobility  of  electrons  in  the  conduction  band  set 
by  lattice  scattering  in  square  centimeters  per  volt-second, 

7.  the  upper  limiting  mobility  of  holes  in  the  valence  band  set  by 
lattice  scattering  in  square  centimeters  per  volt-second, 

8.  the  temperature  in  kelvin  raised  to  the  3/2  power, 

9.  the  relative  effective  mass  of  electrons  in  the  conduction  band, 

10.  the  relative  effective  mass  of  holes  in  the  valence  band, 

11.  the  density  of  states  in  the  conduction  band  in  reciprocal  cubic 
centimeters  (designated  FNC  in  Program  and  on  output  listings), 
and 

12.  the  density  of  states  in  the  valence  band  in  reciprocal  cubic  cen- 
timeters (designated  FNV  in  Program  and  on  output  listings). 

These  parameters  are  printed  across  a  single  line  on  the  output  listing 
in  the  order  given.     Since  these  parameters  are  functions  of  tempera- 
ture,  each  such  line  corresponds  to  a  different  temperature.  The 
resulting  table  is  provided  with  descriptive  column  headings  as  shown  in 
the  sample  output  listing  in  section  5.5.     The  last  five  items  on  the 
above  list  are  useful  for  interpreting  Arrhenius  plots  of  the  logarithm 
of  carrier  density  vs.   1000/T  in  terms  of  the  activation  energies  of  the 
dopant  states  present  (e.g.,   as  might  be  done  in  interpreting  the 
results  of  measurements  of  Hall  coefficient  as  a  function  of  tempera- 
ture ) . 

2.1.4     DOLOOP  XX     ARG(l)  ARG(2) 

This  option  causes  the  program  to  generate  multiple  listings  for  a  sin- 
gle set  of  five  dopant  state  data  cards.     Each  successive  listing  has 
the  density  of  a  selected  dopant  state  incremented  by  a  specified 
amount.     The  following  information  must  be  provided  on  the  DOLOOP  XX 
control  card: 

A.     DOLOOP  XX  is  the  name;  XX  specifies  which  dopant  state  is  to  have 
its  density  incremented.     XX  is  selected  from  the  list  Dl,  D2,  Al, 
A2,   and  A3.     If  anything  else  appears  in  the  XX  position  on  this 
control  card,   an  error  message  to  this  effect  will  be  printed  on  the 
output  listing,   and  the  DOLOOP  XX  control  card  will  be  ignored. 
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B.  ARG(l)  is  a  numerical  argument  which  specifies  the  amount  by  which 
one  would  like  to  increase  (or  decrease)  the  dopant  density  each 
time.     Its  dimensions  are  (centimeters)"^.     In  accordance  with 
the  format  statement  pertaining  to  data  cards  (i.e.,   2A6,  E10.4, 
2F10.4),   it  is  expressed  in  E  format  in  positions  13  to  22  on  the 
control  card. 

C.  ARG(2)  is  a  dimensionless  argviment  which  specifies  the  total  number 
of  listings  desired.     In  accordance  with  the  format  statement 
pertaining  to  data  cards,   it  is  expressed  in  F  format  (with  decimal 
point)  in  positions  23  to  32  on  the  control  card. 

The  DOLOOP  XX  option  refers  only  to  the  next  set  of  five  dopant  state 
data  cards,  and  everything  reverts  back  after  completing  the  number  of 
listings  specified.     Only  one  dopant  state  can  be  incremented  at  a  time 
(i.e.,   no  provision  has  been  made  for  nested  do loops). 

2.1.5  STOP 

This  option  causes  a  FORTRAN  stop  instruction  to  be  executed  and  is  used 
to  exit  the  program.     The  corresponding  control  card  consists  of  the  sin- 
gle word  STOP  in  the  first  four  positions  on  the  card  with  the  rest  of 
the  card  blank.     Normally,   the  STOP  control  card  is  the  last  card  in  the 
deck  of  cards  to  be  read  by  the  program. 

2.2     WARNING  AND  ERROR  MESSAGES 

There  are  several  warning  and  error  messages  that  are  printed  on  the  output 
listing  when  unusual  circumstances  are  encountered  during  the  execution  of 
the  program.     A  warning  message  informs  the  user  that  the  indicated  condition 
was  encountered  and  suggests  that  its  effects  on  the  validity  of  the  next 
data  line(s)  be  considered.     Error  messages  inform  the  user  of  an  event  that 
precludes  the  continuation  of  the  current  run  and,   after  printing  the  error 
message,   the  program  proceeds  to  the  next  set  of  five  dopant  state  data 
cards.     The  messages  are: 

1.  BORN  APPROXIMATION  FOR  IONIZED  SCATTERING  NOT  VALID  FOR  MAJORITY  CARRIER 

This  warning  message  is  printed  when  the  log  term,   BN,   in  the  modified 
Brooks-Herring  ionized  impurity  formula  is  less  than  or  equal  to  10.0. 
In  the  Born  approximation,   the  scattering  is  treated  as  a  small  perturba- 
tion on  the  motion  of  the  incident  carrier.     This  assumption  is  violated 
at  high  carrier  densities  and  low  temperatures. 

2.  (EF-EG)/CAYT  =  XX     ERROR  IN  N  MAY  EXCEED   1  PERCENT 

In  this  warning  message,  XX  is  replaced  by  the  dimensionless  ratio  of  the 
difference  in  energy  between  the  Fermi  level  and  the  conduction  band 
divided  by  the  voltage  equivalent  of  the  temperature  (i.e.,   kT/q) .  This 
message  is  printed  when  XX  >  -2  and  is  a  warning  that  the  calculation  of 
electron  density  may  be  in  error  by  greater  than  1  percent  due  to  the 
expressions  used  to  approximate  the  Fermi -Dirac  integral   [5] .     This  oc- 
curs only  at  high  dopant  densities  when  the  material  is  degenerate  or 
nearly  so. 
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3.  EF/CAYT  =  XX     ERROR  IN  P  MAY  EXCEED   1  PERCENT 

This  warning  message  is  identical  in  function  to  that  of  the  one  immedi- 
ately preceding,   except  that  it  pertains  to  hole  density  instead  of  elec- 
tron density. 

4.  INVALID  ARGUMENT  OF  DOLOOP  COMMAND 

This  error  message  is  printed  when  the  argument  on  the  DOLOOP  XX  control 
card  (i.e.,  XX)  is  not  selected  from  the  list  D1,  D2,  A1,   A2,   or  A3. 
After  printing  this  message,   the  program  proceeds  to  process  the  run 
without  the  DOLOOP  option  (i.e.,   it  produces  one  run  instead  of  the  num- 
ber called  for  on  the  invalid  DOLOOP  XX  control  card). 

5.  CHARGE  BALANCE  EQUATION  COULD  NOT  BE  SOLVED,    IFLAG  =  X 

This  error  message  is  printed  when  the  algorithm  for  the  determination  of 
the  Fermi  level  position  cannot  find  a  solution  to  the  charge  balance 
equation.     When  this  message  is  printed,  X  is  replaced  by  the  value  of  an 
integer  flag  (i.e.,   IFLAG)  that  specifies  the  reason  why  the  charge  bal- 
ance equation  could  not  be  solved.     The  significance  of  the  values  of 
IFLAG  are  discussed  in  section  3.1  below. 

6.  LOOPSWITCH  IS  OUT  OF  RANGE 

This  is  an  error  message  that  should  never  be  printed  if  the  program  is 
functioning  properly.     LOOPSWITCH  is  a  program  variable  whose  value  is 
determined  by  the  XX  portion  of  the  DOLOOP  XX  control  card  and  the  pro- 
gram does  not  contain  a  statement  capable  of  causing  it  to  be  out  of 
range.     This  message  was  included  as  a  check  on  the  proper  operation  of 
the  DOLOOP  option. 

2.3     ADDITIONAL  INFORMATION 

The  program  listing  in  section  4  of  this  report  is  copiously  documented  with 
comments  to  guide  the  user  through  the  details  of  the  program.     The  complete 
program,   listed  in  section  4,   consists  of  a  main  program  followed  by  a  series 
of  three  subroutines  and  one  function.     The  function  TEST  and  the  subroutine 
ZEROIN  are  used  for  finding  the  position  of  the  Fermi  level,   and  the  subrou- 
tines MOB  and  SICI  are  used  for  performing  the  mobility  computations.  Notice 
that  the  mobility  algorithm  requires  knowledge  of  many  of  the  parameters  cal- 
culated by  the  Fermi  level  locating  routine;  thus,   it  cannot  be  used  without 
first  finding  the  location  of  the  Fermi  level. 

All  of  these  routines  have  been  written  in  FORTRAN  IV  for  use  on  the  Univac 
1108  computer  at  the  National  Bureau  of  Standards.     This  is  a  36-bit  machine 
with  26-bit  mantissa  floating-point  numbers.     It  may  be  necessary  to  use 
double  precision  to  obtain  a  satisfactory  accuracy  on  machines  with  signifi- 
cantly smaller  word  lengths.     It  may  also  be  necessary  to  modify  the  alpha- 
numeric format  for  data  cards  and  the  IF  statements  that  test  alphanumeric 
variables  for  machines  that  do  not  store  six  alphanumeric  characters  per 
word. 
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There  are  several  places  within  the  program  where  a  zero  denominator  might 
occur  during  execution  (e.g.,  when  computing  the  Hall  coefficient  for  a  very 
small  carrier  density  at  low  temperatures).     The  present  program  tests  for 
this  circumstance  in  all  places  where  this  has  occurred  in  actual  use.  When 
this  condition  is  encountered,   a  suitable  flag  is  set,   the  division  in  ques- 
tion is  circumvented,   and  the  printing  of  all  affected  parameters  is  sup- 
pressed.    In  addition,   all  parameters  with  a  zero  value  are  omitted  from  the 
output  listing  so  as  not  to  clutter  the  output  pages  with  multiple  zero  en- 
tries.    This  has  been  accomplished  by  preceding  each  WRITE  statement  with  an 
appropriate  IF  statement  that  suppresses  the  printing  of  undesired  output. 
Since  this  may  slow  down  the  printing  speed  or  increase  the  cost  of  execution 
in  some  systems,   the  reader  may  wish  to  modify  this  portion  of  the  program  to 
achieve  the  same  result  with  object-time  format  statements   [5]  computed 
during  execution  and  stored  in  an  array  for  use  by  the  appropriate  WRITE 
statement. 

Copies  of  the  main  program  are  available  on  cards  from  the  authors. 

3.     THEORY  AND  METHOD  OF  COMPUTATION 

The  computations  of  the  position  of  the  Fermi  level  and  carrier  mobilities 
are  key  parts  of  this  program.     The  various  equations  used  in  these  computa- 
tions are  discussed  below.     Numbers  cited  in  the  following  sections  refer  to 
statement  numbers  in  the  complete  program  given  in  section  4.     If  two  numbers 
are  given,   the  first  refers  to  electrons  or  donor  states,   and  the  second  to 
holes  or  acceptor  states. 

3.1     LOCATION  OF  THE  FERMI  LEVEL 

For  uniformly  doped  silicon  in  thermal  equilibrium  (as  assumed  herein),  the 
summation  of  all  negative  charges  must  equal  the  summation  of  all  positive 
charges  (i.e.,   the  material  must  be  electrically  neutral).     The  negative 
charges  consist  of  electrons  in  the  conduction  band  and  any  negatively 
charged  impurity  states  (acceptors).     The  positive  charges  consist  of  holes 
in  the  valence  band  and  any  positively  charged  impurity  states  (donors).  The 
Fermi  level  is  located  at  just  that  point  that  makes  the  total  net  charge 
equal  zero.     The  FORTRAN  parameter  TEST  defined  in  statement  140  of  the  func- 
tion TEST  is  a  measure  of  how  precisely  the  algorithm  has  located  the  posi- 
tion of  the  Fermi  level  to  obtain  charge  neutrality. 

TEST  =  (Net  negative  charge)  -  (Net  positive  charge) 

TEST  =   (CN  +  DI3  +  DI4  +  DI5)   -   (CP  +  DI 1   +  DI2)  (140)* 

where  the  symbols  represent  densities  in  reciprocal  cubic  centimeters  as  fol- 
lows : 

CN,  electrons  in  the  conduction  band, 

CP,  holes  in  the  valence  band, 

DI 1 ,  ionized  donor  species  1, 


*The  statement  numbers  in  this  section  refer  to  the  function  TEST  unless 
otherwise  indicated. 
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DI2,  ionized  donor  species  2, 

DI3,  ionized  acceptor  species  3, 

DI4,  ionized  acceptor  species  4,  and 

DI5,  ionized  acceptor  species  5. 

The  densities  of  electrons  and  holes  in  the  conduction  and  valence  bands, 
respectively,  are  computed  using  Blakemore's  approximations  to  the  Fermi- 
Dirac  integral   [5] : 

n  =  N   [exp(-  n)  +  0.27]"^  (8,26) 
for  values  of  n  in  the  range  from  -80  to  +1, 

2  3/4 

n  =  0.752253  N  (n    +  1.7)^  (10,30) 
for  values  of  n  greater  than  +1,  and 

n  =  0  (2,20) 

for  values  of  ri  less  than  -80.     In  these  equations,  n  is  the  electron  or  hole 
density,  N  is  the  density  of  states  in  the  appropriate  band,    n  =  q(Ep  - 
Eg) At  for  electrons  (4),  and  -qEp/kT  for  holes  (22),  Ep  is  the  Fermi 
energy,   in  electron  volts.  Eg  is  the  width  of  the  forbidden  energy  gap  in 
electron  volts  (the  zero  of  energy  is  assumed  at  the  valence  band  edge),  and 
kT/q  is  the  voltage  equivalent  of  the  temperature  (353  of  Main  Program). 

The  density  of  states  is  given  by 

2  3/2     *3/2  3/2 
N  =  2(2tt  k  m  A  )         m  1 
o 


A  ooo  *3/2  3/2 

=  4.829  X  10       m  T 


(354,355)  Main  Program 


where  m    is  the  effective  mass  of  the  appropriate  carrier  relative  to  the 
free  electron  mass,  m^.     The  relative  electron  effective  mass,  m^,  is 
given  by  [7]''" 

m*  =   1.0627  -   1  .61708  x  lO'^^T  +  6.83008  x  10~^T^  -  3.32013  x  10"^"^ 
e 

-11  4  -14  5  -17  6 

+  8.04032  X  10       T     -  9.66067  x  10       T     +  4.54649  x  10  T 

(350)  Main  Program 

and  the  relative  hole  effective  mass,  m^,   is  given  by  [7] 


This  and  the  following  two  equations  have  been  carried  to  high  order  in 
kelvin  temperature  for  purposes  of  computational  accuracy,  without  regard  to 
how  high  an  order  is  really  necessary  to  describe  the  measured  data  to 
within  experimental  accuracy. 
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m*  =  0.590525  -  5.23548  x  ^0~\  +  1.85678  x  10~^T^  -  9.67212  x  10~\"^ 
h 

+  2.30049  X  10~^°T^  -  2.596730  x  10~^\^  +  1.11997  x  lo'^^T^  . 

(352)  Main  Program 

The  temperature  dependence  of  the  forbidden  energy  gap  of  silicon  is  taken 
into  account  in  accordance  with  a  fit  to  the  data  of  Macfarlane  et  at.   [8] : 

-5  -7  2  -10  3 

E     =  1.15556  +  3.23741  x  10     T  -  8.7011  x  10     T     +  9.95401  x  10  T 

g 


-  3.80977  X  10~^^T^ 


(356)  Main  Program 


The  density  of  ionized  impurities  is  computed  from  the  following: 

N.  =  N     [1  +  g  exp(A)]  (40  ff) 

It 

where        is  the  total  density  of  the  dopant  impurity,   g  is  the  degeneracy 
factor,   and  A  =  q(Ep  -  Eg  +  E^)AT  for  donor  states  and  q(E^  -  Ep)AT 
for  acceptor  states,  and  E^^  is  the  value  of  the  activation  energy,  relative 
to  the  appropriate  band  edge.     This  treatment  of  ionized  impurity  states  is 
only  valid  for  simple  donor  and  acceptor  states  that  have  the  following  prop- 
erties: 

1.  Each  impurity  state  gives  rise  to  only  one  electronic  state.  The 
present  version  of  this  program  cannot  handle  the  case  of  coupled 
multiple  levels  associated  with  a  single  impurity  species  unless  the 
levels  are  widely  separated  and  can  be  treated  independently. 

2.  The  donor  levels  are  neutral  when  occupied  by  an  electron  and  have  a 
temperature-independent  activation  energy  referenced  to  the  bottom 
of  the  conduction  band. 


3.  The  acceptor  levels  are  neutral  when  occupied  by  a  hole  and  have  a 
temperature-independent  activation  energy  referenced  to  the  top  of 
the  valence  band. 

These  properties  are  appropriate  for  the  colximn  3  and  column  5  impurities  in 
the  periodic  table  (e.g.,  boron  and  phosphorus)  and  may,   or  may  not,  be 
appropriate  for  other  impurities  with  larger  activation  energies  or  multiple 
energy  levels.     Impurity  states  satisfying  these  criteria  are  called  dopant 
states  in  this  report  and  are  characterized  by  the  following  three  parameters 
which  are  specified  on  the  input  data  cards  for  each  state: 

1.  the  concentration  of  the  dopant  state  per  unit  volume  of  silicon, 

2.  the  limiting  activation  energy  at  low  concentrations  relative  to  the 
appropriate  band  edge,  and 

3.  the  degeneracy  of  the  dopant  state  (an  integer  larger  than  unity). 

The  decrease  in  activation  energy  with  increasing  dopant  density  or  increas- 
ing compensation  is  taken  into  account  with  the  use  of  an  equation  based  on 
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the  work  of  Penin  et  at.    [9]   for  phosphorus-doped  silicon.     Penin's  expres- 
sion gives  zero  activation  energy  at  3  x  10^^  cm~-^  regardless  of  the 
initial  activation  energy,   E^g •     Since  deeper  levels  are  expected  to  re- 
quire greater  densities  to  obtain  zero  activation  energy  [10],  Penin's  ex- 
pression has  been  modified  to  retain  the  same  decrease  in  activation  energy, 
in  eV,  for  a  given  increase  in  dopant  density.     To  accomplish  this,  Penin's 
factor  of  8  X  0.045,  where  0.045  is  the  phosphorus  activation  energy,  is 
fixed  at  0.36.     Thus,   since  E^^q        greater  for  a  deep  level,   the  dopant 
density  for  zero  activation  energy  is  correspondingly  greater. 

The  resulting  equation  for  calculating  the  decrease  in  activation  energy  with 
dopant  density  and  compensation  is 

=  -  +  7")  (-r  ,/a*)  -  a  , 

A         AO  2       r  d  2 

d 

(70  ff)  Main  Program 

where  E^q  is  the  activation  energy  for  low  densities  of  the  dopant,  = 
[3/{47rNi  )  ]         is  the  average  distance  between  the  dominant  dopant  atoms  of  a 
given  class  (donor  or  acceptor),   a*  =  21  x  10"^  cm  is  the  effective  Bohr 
radius,  a  =  3.6  x  10"^  eV  cm,   and  N2  is  the  total  density  of  compensating 
impurities  of  the  other  class  (acceptor  or  donor).     This  ec[uation  is  used  in 
the  computer  program  for  all  dopant  states,  but  its  applicability  to  deep 
levels  has  not  been  verified. 

Since  the  reduction  in  activation  energy  with  increasing  donor  density  is  due 
in  large  part  to  the  lowering  of  the  bottom  of  the  conduction  band,   the  total 
donor  density  is  used  in  the  calculation  of  the  reduction  for  each  donor 
state.     Likewise,  when  more  than  one  acceptor  state  is  present,   the  total 
acceptor  state  density  is  used  to  calculate  the  reduction  in  activation  ener- 
gy of  each  acceptor  state.     However,  when  two  or  more  donors  (or  acceptors) 
are  present,   each  in  significant  quantity,   this  approach  probably  overesti- 
mates the  reductions  in  activation  energy  of  each  acceptor  state. 

The  computation  of  the  reduction  in  activation  energy  is  not  included  as  part 
of  the  routine  for  finding  the  Fermi  level,  but  has  been  included  in  the  Main 
Program  as  part  of  the  initialization  procedure.  The  routine  for  finding  the 
Fermi  level  position  assumes  that  the  values  of  activation  energy  previously 
placed  in  COMMON  are  the  correct  values.  The  complete  program  outputs  both 
the  limiting  value  of  activation  energy  initially  inputted  and  the  calculated 
lowered  value  of  activation  energy  (see  sec.  5). 

The  correct  position  of  the  Fermi  level  is  found  by  the  subroutine  ZEROIN 
[11-13]  which  searches  for  a  zero  of  the  parameter  TEST.     ZEROIN  is  listed  as 
part  of  the  complete  program  in  section  4.     The  subroutine  is  well  documented 
with  descriptive  comments  before  the  code  and  explanations  of  all  significant 
steps  as  they  occur  in  the  body  of  the  code.     The  subroutine  looks  for  a  zero 
between  the  variables  EF  and  EFG  which  are  given  initial  values  of  0.0  and 
1.2  eV,   respectively,  by  the  calling  program.     In  searching  for  the  root, 
ZEROIN  uses  the  secant  rule  unless  tests  indicate  that  bisection  would  be 
advantageous.     The  interval  between  EF  and  EFG  is  narrowed  until  the  stopping 
criterion  is  satisfied.     The  stopping  criterion  is 
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EF-EFG 


<  RE  X  EF  +  AE  , 


where  RE  is  the  maximum  allowable  relative  error  and  AE  is  the  maximum  allow- 
able absolute  error  in  eV.     In  the  calling  program  these  are  both  assigned 
the  value  1  x  10"^  so  that  they  have  equal  weight  in  the  stopping  crite- 
rion when  EF  =  1.0  eV.     Initially,   the  subroutine  compares  RE  with  the  round- 
off error  of  the  machine,  ER,   specified  in  a  DATA  statement,  which  can  be 
easily  changed  to  correspond  to  the  user's  computer,   and  selects  the  larger 
of  the  two  for  the  stopping  criterion. 


At  the  conclusion  of  the  computations,  EF  is  the  better  approximation  to  the 
zero  of  TEST  as  the  subroutine,   by  interchange  if  necessary,  makes 
I  TEST(EF)  I  <  I  TEST(EFG)  |  .     Before  the  return  to  the  main  program,   the  param- 
eter IFLAG  is  set  to  indicate  the  status  of  the  results.     The  normal  case 
with  the  root  located  within  the  stopping  criterion  is  denoted  by  IFLAG  =  1. 
A  possible,  but  unlikely,   result  is  TEST(EF)  =  0  without  the  stopping  crite- 
rion being  met  (current  value  of  EFG  not  close  enough  to  EF ) .     For  this  situ- 
ation IFLAG  =2.     If  the  stopping  criterion  is  met,  but  |  TEST(EF)  |  is  larger 
than  the  absolute  value  of  TEST  evaluated  with  either  initial  argument,  IFLAG 
=3.     If  the  stopping  criterion  is  satisfied,   but  TEST(EF)  x  TEST (EFG)   >  0, 
there  is  apparently  no  root  in  the  interval  and  IFLAG  =  4.     ZERQIN  permits  up 
to  500  evaluations  of  TEST  before  terminating  and  setting  IFLAG  =5.     If  upon 
return  to  the  main  program  IFLAG  is  not  1  or  2,   its  value  is  printed  for  di- 
agnostic purposes.     The  program  then  proceeds  to  the  next  temperature.     For  a 
successful  evaluation  (IFLAG  =  1  or  2),   the  residual  value  of  TEST,  called 
TESTF,   is  printed  in  the  DOUBLE  ON  output  listing. 

3.2     COMPUTATION  OF  CARRIER  MOBILITIES  AND   RELATED  PARAMETERS 


The  scattering  mechanisms  which  are  included  in  the  calculation  of  the  carri- 
er mobility  are  lattice,   ionized  impurity,   neutral  impurity,   and  electron- 
electron  (or  hole-hole)  interactions.     The  mobility  calculations  are  sum- 
marized in  the  tables  which  follow.     In  table  1,   the  symbols  used  in  the 
mobility  formulas  are  defined  and  the  corresponding  computer  program  name  is 
given.     Table  2  gives  the  mathematical  formulas  for  calculating  electron 
mobility  along  with  the  statement  number  of  the  corresponding  expression  in 
the  computer  program.     Table  3  gives  this  information  for  the  calculation  of 
hole  mobility. 

The  temperature  dependence  of  the  electron  mobility  due  to  lattice  scattering 
is  taken  from  the  work  of  Norton  et  al.    [14],   and  the  magnitude  is  based  on 
the  results  at  300  K  of  Thurber  et  at.    [15].     The  temperature  dependence  of 
this  mobility  changes  slope  at  about  110  K,   and  this  is  the  reason  for  the 
different  equations  in  table  2.     For  holes,   the  lattice  mobility  for  tempera- 
tures below  72  K  is  from  Braggins   [15],  whereas  at  higher  temperatures  the 
results  of  Ludwig  and  Watters  [17]  are  used. 

The  ionized  impurity  scattering  mobility  is  calculated  by  the  Brooks-Herring 
[18,19]  expression,  with  the  modifications  discussed  by  Li  and  Thurber  [20] 
and  Li   [21],     For  p-type  material,   individual  mobilities  are  computed  for 
each  band,   and  then  the  total  mobility  is  found  by  weighting  the  mobility  in 
each  band  by  the  hole  density  in  that  band.     The  effect  of  carrier-carrier 
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Table  1. 


Symbols  Used  in  Mobility  Formulas. 


rid  L itcllla  L  l  Ud  i 

Symbol 

r  1  u  y  1  uM  t 

Name 

Definition* 

u 

SUMTD 

Total  donor  density  (cm" 3) 

SUMTA 

Total  acceptor  density  (cm"^) 

III  N 

U  L 11 

I  ^ff'\r0  mnhil'il'v/  "Frip  plpp1">"nn<;    ( rvn^  /\!  • 

ULP 

Lattice  mobility  for  holes  (cm^/V-s) 

UIN 

Ionized  impurity  mobility  for  electrons  (cm^/V-s) 

1 1  T  P 

Tnm*  7pH   imniir^'i       mnh"i  1  i  I'v  "for  h^lp<^  irm^/V*'^) 

Y 

GAMMA 

Factor  giving  effect  of  carrier-carrier  scattering 
on  ionized  impurity  mobility 

1 

CION 

Density  of  ionized  impurities  (cm" 3) 

G(b) 

GN 

Log  term  in  Brooks-Herring  ionized  impurity  mobility 
formula  for  electrons 

U  r 

Inn   "hp  v^m   "in   Rv^nnl^Q  —  Mpv^K'inn    inni'yprl   TmniiP"it"V/  mn  hi  1  i  "I"  v 

formula  for  holes 

SUMID 

Total  ionized  donor  density  (cm" 2) 

V 

SUMIA 

Total  ionized  acceptor  density  (cm" 3) 

n' 

PN^TAR 

LiliO  1  r\r\ 

F*F"f pr "t"  1  \/p  plprtK^nn  <^r^ppninn  Hpn*^i1"\/  frm"^^ 

L.  \   ICl-LIVC     CICV^LlUll     OV^ICCIIIIiy     UCIIOIUjf      ^V-IH  J 

P 

C"P*fQ^"t""iv/£i    ki/~\np    c^v^ppnnnn  HpnCT'hx/ 

tTTcCLiVc  no  1  c  bcrcciiiiiy  uciibiLy  ^lim  j 

n  T  Fl 

tTT  LI  mes  relative  qiciccliil  LuiisLaiiu  Lniicb  \jc  i  - 
mittivity  of  free  space 

k 

Boltzmann  constant 

h 

Planck  constant  (li  =  h/2-n) 

q 

Electronic  charge 

mo 

Free  electron  mass 

< 

Longitudinal  effective  mass  (in  units  of  mp) 

* 

Densi ty-of-states  effective  mass  (in  units  of  m^) 

* 

Conductivity  effective  mass  (in  units  of  m^) 

UNIS 

Neutral  impurity  scattering  mobility  (cm^/V-s) 

\ 

CNUT 

Density  of  neutral  impurities  (cm"^) 

n 

CN 

Electron  density  (cm"^) 

P 

CP 

Hole  density  (cm" 3) 

*  Quantities  in  SI  units  unless  otherwise  stated. 
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Table  2. 


Formulas  for  Calculating  Electron  Mobility. 


Scatteri  ng 
Mechani  sm 

Formula 

Sta  tement 
Number 

1 fl  f ti  TP  + 

electron-electron 

Nn  <  3  X  IQl^  cm-3 

111]         ^          ij                          1    \J  Lrlll 

T  <  108  K    ui   =  1  89  X  107T-1-52  ri  _  O  08  Nn/f2  x 

1017)1 
1  u     ;  J 

Q 
O 

T  >  108  K    M|_  =  7.98  X  108T-2-32  [i  _  o.08  Nd/(2  x 

1017)] 

10 

Nd  >  3  X  1017  cm-3 

T  £  108  K    ML  =  l -89  X  107T-1-52  (0.884) 

22 

T  >  108  K    vi|_  =  7.98  x  io8t-2-32  (0.884) 

30 

Ionized  impurity  + 
el ectron-electron 

yj  =  Y  X  7.3  X  lOl^jl. 5/NiG(b) 

Y  =        [1  -  exp(-n'/ND+)] 

n'  =  n  +  p  +  Nd+  (1  -  Nd+/Nd) 

G(b)  =  ln(b  +  1 )  -  b/(b  +  1 ) 

24k  e    m*  (kT)2 
b  =          °    ^             X  10-6 
q2li2n' 

90 

86 

42 
80 

52 

Combined  lattice 
and  ionized  im- 
purity 

P|  T    =   Pi  f(X) 
LI  L 

f(X)  =  1  +  X2  [CiXcosX  +  sinX(SiX  -  J)] 
X2  =  6 

104 
104 
98 

Neutral  impurity 

_  n  no         r2   fkTvl/2       1    , Em  1/2 

P    _  1.136  X  10-19  m^ 
t..     —       —  ...... 

114 
112 

k2  mo 

_     2TT3n3  mi 

no 

Total  electron  y  =  (Vp,  t  +  l/y.,)"i  116 

mobility  ^ 


*These  statement  numbers  refer  to  subroutine  MOB. 
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Table  3. 


Formulas  for  Calculating  Hole  Mobility. 


Scattering 
Mechanism 


Formula 


Statement 
Number* 


Lattice  + 

N„  <  3  X  1017  cm-3 
A 

hole-hole 

T  < 

72  K        =  1.6  X  lo^T-i-S^  [1  -  0.08  Na/(2  x  lO^^)] 

118 

T  > 

72  K        =  2.3  X  lon-2-7  [1  .  0.08  Na/(2  x  lO^^)] 
>_  3  X  1017 

120 

T  < 

72  K        =  1.6  X  107T-1-54  (0.884) 

132 

1  > 

79    V             -    9    "3         in9T-2    7           DQ/I  ^ 
1 C   K           ~      • lU    1      •  ^U.cScS'tj 

1  /I  n 
1  f  u 

Ionized  impurity 

+  hnl p-hnl p 

yCpj^  +        ('"d2/'"dl  ^ ^ ■ '     ^13  (md3/>)^-^exp(-0.044/kT)] 

1 96 

[1  +  (mj3/mJ^)i-5  +  (m33/m;^)i-5  exp(-0.044/kT)] 
[1  -  exp(-pVNj)] 

Y  = 

198 

P'  = 

p  +  n  +  N+  (1  -  N+/N„) 
^              A            A  A 

152 

for 

i  =  1,  2,  3 

^li 

23-5  (4TTKe  )2  m*.0-5(kT)i-5 

0         dl                        ^  ^Q_2 

ui-5  q3  m*.  Nj  G(b.) 

192 

G(b. 

)  =  ln(b.  +  1)  -  b./(b.  +  1) 

190 

24Ke    m*.  (kT)2 
0  di 

q2ti2p' 


10-5 


162 


Combined  lattice 
and  ionized  im- 
puri  ty 


-  f(X) 

f(X)  =  1  +  X2[CiXcosX  +  sinX(SiX  -  J)] 


X2  =  6  P^^^I 


214 

214 
208 


Neutral  impurity 


,(m*  /m*  )i-5  +  exp(-0.044/kT)] 


Nl      ^N2^"'d2""dl^       '  ^N3^"'d3""dl 


^N 


[1  +  (m*2/mdi)^-5  +  (md3/'^dl)^"^  exp(-0.044/kT)] 


for  i  =  1 ,  2,  3 
^Ni  =  [|  (|I-) 


1/2 


.     E  1/2 

+  -  (^) 

3  ^  kT^ 


1.136  X  10-19  m*. 


<^  m. 


2Tr3q  3m 


di 


^"^      5N.,4TrKe  h3 
N  0 


X  10- 


224 


222 


221 


220 


Total  hole 
mobil  ity 


p  =  (I/pli  +  Vp,^)-i 


226 


These  statement  numbers  refer  to  subroutine  MOB. 
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scattering  in  reducing  the  lattice  and  ionized  impurity  mobilities  is  taken 
into  account  following  the  work  of  Li  and  Thurber   [20].     The  mobility  due  to 
both  lattice  and  ionized  impurity  scattering  is  combined  according  to  the 
mixed-scattering  formula  of  Debye  and  Conwell   [22]   involving  sine  and  cosine 
integrals. 

The  expressions  used  to  calculate  the  carrier  mobilities  contain  the  dielec- 
tric constant  of  silicon  as  a  parameter.     The  value  of  11.7  for  the  relative 
dielectric  constant  of  silicon  at  300  K  was  taken  from  Dunlap  and  Watters 
[23].     The  temperature  variation  was  calculated  from  the  index  of  refraction 
data  of  Cardona  et  at.    [24],   assuming  that  the  index  of  refraction  is  the 
square  of  the  relative  dielectric  constant.     The  resulting  expression  for  the 
program  parameter  DIEL  evaluated  in  subroutine  MOB  statement  No.  2  is: 

DIEL  =  4     (absolute  dielectric  constant)  and 

DIEL  =  (1.2711  X  10"^)  exp(7.8  x  10~^T)  farads  per  meter. 

The  neutral  impurity  scattering  mobility  is  calculated  from  the  Erginsoy 
equation  [25] ,  with  the  modification  proposed  by  Sclar  [26] .     The  net  mobil- 
ity is  found  by  combining  reciprocally  the  neutral  impurity  mobility  with  the 
result  for  the  combined  lattice  and  ionized  impurity  mobility. 

Both  electron  and  hole  mobilities  are  calculated,   even  though  usually  only 
one  carrier  is  important.     An  effective  conductivity  mobility  is  calculated 
by  summing  the  electron  and  hole  mobilities  weighted  by  the  density  of  each 
carrier.     The  effective  conductivity  mobility  computed  this  way  is  the  same 
as  the  majority  carrier  conductivity  mobility  in  extrinsic  specimens.  In 
intrinsic  specimens,   it  corresponds  to  a  hypothetical  mobility  for  the  total 
carrier  density  (i.e.,   electron  plus  hole  density).     The  effective  conductiv- 
ity mobility  is  the  only  value  of  mobility  listed  in  the  no-option  output 
format.     Additional  mobility  values  are  available  in  the  DOUBLE  ON  output 
format.     The  expressions  for  calculating  the  conductivity,   resistivity,  and 
Hall  coefficient  include  appropriately  weighted  contributions  from  both 
charge  carriers. 
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4.     Annotated  Program  Listing 


C  ELECTRICAL  PROPERTIES  OF  SILICON 
C 

C  PROGRAM  CAN  HANDLE  2  DONOR  LEVELS  AND  3  ACCEPTOR  LEVELS 

C 

C 

C  VARIABLE  NAMES 
C 

C  LEVEL  NAME  DENSITY        ENERGY  DEGEN.  ION.  CONC. 

C 

C  DONOR  1         NAME1  DNI  ENINI  DGl  DI1 

C  DONOR  2         NAME2  DN2  ENIN2  DG2  DI2 

C  ACCEPTOR  1  NAME3  DN3  ENIN3  DG3  DI3 

C  ACCEPTOR  2  NAME4  DN4  ENIN4  DG4  DI4 

C  ACCEPTOR  3  NAMES  DN5  ENIN5  DG5  DI5 

C 
C 

C  STATEMENT  NUMBERS  ENDING  IN  ZERO  ARE  PART  OF  THE  PROGRAM  LOGIC 
C  OTHER  NUMBERS  WERE  ADDED  FOR  IDENTIFICATION  PURPOSES 

C  A  MODIFIED  MKS  SYSTEM  OF  UNITS  IS  USED  WITH  CM  REPLACING  METERS 
C 

C  T  =  THE  TEMPERATURE  IN  KELVIN 

C  EG  =  THE  ENERGY  GAP  OF  SILICON  IN  ELECTRON  VOLTS 

C  EXCITON  ENERGY  OF  0.01   EV  IS  NOT  INCLUDED  IN  EG 

C  EF  =  FERMI  LEVEL  IN  EV  (DETERMINED  IN  PROGRAM  BY  TRIAL  AND  ERROR) 

C  EF  IS  TAKEN  TO  BE  ZERO  AT  THE  VALENCE  BAND  EDGE 

C  CN  =  THE  CONCENTRATION  OF  ELECTRONS  IN  THE  CONDUCTION  BAND  IN  CM-3 

C  CP  =  THE  CONCENTRATION  OF  HOLES  IN  THE  VALENCE  BAND  IN  CM-3 

C  MOB  =  MOBILITY  ON  THE  OUTPUT  PAGES,  BUT  NOT  IN  THE  PROGRAM: 

C  IN  THE  PROGRAM,  THE  FOLLOWING  NOTATION  IS  USED: 

C  UDN  =  THE  ELECTRON  MOBILITY  IN  CM**2/ (VOLT-SEC ) 

C  UDP  =  THE  HOLE  MOBILITY  IN  CM**2/ (VOLT-SEC) 

C  UD  =  THE  CONDUCTIVITY  MOBILITY  IN  CM**2/ (VOLT-SEC ) 

C  RH  =  THE  HALL  COEFFICIENT  IN  CM**3/C0UL0MB 

C  SIGMA  =  THE  CONDUCTIVITY  IN  MHO/CM 

C  RHO  =  THE  RESISTIVITY  IN  OHM-CM 

C  ON  THE  OUTPUT  PAGE  THE  FOLLOWING  NOTATION  HAS  BEEN  USED: 
C  Dl   ION  =  THE  DENSITY  OF  IONIZED  DONOR  1   IN  CM-3 

C  D2  ION  =  THE  DENSITY  OF  IONIZED  DONOR  2  IN  CM-3 

C  A1   ION  =  THE  DENSITY  OF  IONIZED  ACCEPTOR  1   IN  CM-3 

C  A2  ION  =  THE  DENSITY  OF  IONIZED  ACCEPTOR  2  IN  CM-3 

C  A3  ION  =  THE  DENSITY  OF  IONIZED  ACCEPTOR  3  IN  CM-3 

C  FNC  AND  FNV  ARE  THE  DENSITY  OF  STATES  IN  THE  CONDUCTION  AND 
C  VALENCE  BANDS  RESPECTIVELY  IN  CM-3 

C 

C  INPUT  DATA  CARDS: 

C  ORDER  OF  CARDS  IS  THE  SAME  AS  IN  THE  TABLE  ABOVE 

C  ORDER  OF  DATA  ON  CARD  =  NAME,  DENSITY,  ENERGY,  DEGEN. 

C  TO  ELIMINATE  A  LEVEL  FROM  THE  PROGRAM,  NAME  IT  'DUMMY' 

C  THE  LAST  DATA  CARD  SHOULD  READ  SIMPLY  'STOP' 

C 

C  THE  FIRST  DATA  CARD  OF  ANY  SET  CAN  BE  USED  TO  CONTROL  THE  PROGRAM 

C  SEE  THE  LIST  OF  IF  STATEMENTS  FOLLOWING  STATEMENT  10  FOR  SUCH  USE 
C  'DOUBLE  ON'  TURNS  ON  2  OUTPUT  DATA  LINES  AT  EACH  TEMP. 
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C  'DOUBLE  OFF'  RESTORES  SINGLE  OUTPUT  LINE  FORMAT 

C  'CONSTANTS'  GENERATES  A  PAGE  OF  THE  PARAMETERS  THAT  DO  NOT 

C  CHANGE  WITH  DOPING  LEVEL 

C  'DOLOOP  XX'    (WHERE  XX  =  D1   -  A3)  PRODUCES  MULTIPLE  RUNS 
C  CARD  FORMAT  =  DOLOOP  XX,  INC.   DENSITY,  NO.  OF  RUNS 

C  'STOP'   IS  USED  TO  TERMINATE  THE  SESSION  WITH  THE  COMPUTER 
C 


EXTERNAL  TEST 

DIMENSION  NAME1 (2) ,NAME2 (2) ,NAME3(2) ,NAME4(2) ,NAME5 (2) 

COMMON  NAME1 ,NAME2,NAME3,NAME4, NAMES 

COMMON  DN1 ,DN2,DN3,DN4,DN5,SUMTD,SUMTA 

COMMON  DI1,DI2,DI3,DI4,DI5 

COMMON  EN1,EN2,EN3,EN4,EN5 

COMMON  DG1,DG2,DG3,DG4,DG5 

COMMON  T,CN, CP, EG, EFME,EFMH,FNC,FNV,CAYT,TPOWR,Q, ALPHA, BETA 
COMMON  UD,UDN,ULN,UDP,ULP,SIGMA,RHO,RH 
COMMON  0UT8,0UT11,ISTAT 
Q=1.6022E-19 

C  INITIALIZE  THE  SETTINGS  OF  THE  OPTION  SWITCHES 

LOOPSW=0. 
D0UBSW=6H0FF 

C 

C  READ  THE  Dl   DATA  (OR  CONTROL)  CARD 

10        READ  20,  NAME1,DN1,ENIN1,DG1 

IF  (NAME1 (1) .EQ.6HD0UBLE)  GO  TO  730 

IF  (NAME1 (1) .EQ.6HC0NSTA)  GO  TO  50 

IF  (NAMEKI)  .EQ.6HD0L00P)  GO  TO  740 

IF  (NAME1 (1).EQ.6HST0P     )  GO  TO  760 

C 

C  READ  THE  REMAINING  FOUR  DATA  CARDS 

READ  20,  NAME2,DN2,ENIN2,DG2 

READ  20,  NAME3,DN3,ENIN3,DG3 

READ  20,  NAME4,DN4,ENIN4,DG4 

READ  20,  NAMES, DNS, ENIN5,DG5 
20        FORMAT  (2A6,E10.4,2F10.4) 
C 

IF  (LOOPSW.EQ.O.)  GO  TO  50 
C  PERFORM  THE  NECESSARY  DOLOOP  INSTRUCTIONS 

DO  720  NRUN=1,NDIT,1 
IF  (NRUN.EQ.1.)  GO  TO  SO 
IF  (L00PSW.EQ.1.)  DN1=DN1+DELDEN 
IF  (L00PSW.EQ.2. )  DN2=DN2+DELDEN 
IF  (L00PSW.EQ.3.)  DN3=DN3+DELDEN  ^ 
IF  (L00PSW.EQ.4.)  DN4=DN4+DELDEN 
IF  (LOOPSW.EQ.S.)  DN5=DNS+DELDEN 
IF  (LOOPSW.GT.S.)  GO  TO  30 
GO  TO  SO 


30        PRINT  40 

40        FORMAT  (1 H0,26HLOOPSWITCH  IS  OUT  OF  RANGE) 
GO  TO  10 

C  RUNN0=1.   DURING  FIRST  RUN,  THEN  2.  UP  TO  77  KELVIN,  THEN  3.  TO 

C  ROOM  TEMPERATURE,  THEN  4.  ABOVE  ROOM  TEMPERATURE 

50        RUNN0=1 . 
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C  CALULATE  REDUCTION  IN  ACTIVATION  ENERGIES 

C  EQUATION  OF  RENIN  ET  AL.   FOR  PHOSPHORUS-DOPED  SILICON  IS 

C  GENERALIZED  FOR  USE  WITH  OTHER  ACTIVATION  ENERGIES 

C  SEE  PAPER  IN  SOVIET  PHYS. -SOLID  STATE  7,  2580  (1966) 

SUMTD=DN1+DN2 

C  GUARD  AGAINST  ZERO  DIVISOR  IN  CALCULATION  OF  RDD 

IF  (SUMTD.LE.O.)  RDD=1.E20 

IF  (SUMTD.LE.O.)  GO  TO  60 

RDD=1 .E8*(3./(4.*3.1416*SUMTD))**0.3333 
60  SUMTA=DN3+DN4+DN5 

C  GUARD  AGAINST  ZERO  DIVISOR  IN  CALCULATION  OF  RDA 

IF  (SUMTA.LE.O.)  RDA=1.E20 

IF  (SUMTA.LE.O.)  GO  TO  70 

RDA=1 . E8* (3. / (4. *3 . 1 41 6*SUMTA) )**0. 3333 
70  EN1=ENIN1 

C  FOLLOWING  STATEMENT  COMPUTES  THE  REDUCTION  IN  ACTIVATION  ENERGY 

C  OMIT  THIS  STATEMENT  IF  FIXED  ENERGY  LEVEL  ENIN1   IS  DESIRED 

C  SAME  COMMENTS  APPLY  TO  OTHER  DOPANT  STATES 

EN1 =ENIN1 -0.36* (0.5+21 ./RDD )*EXP(-RDD/21 .)-3.6E-8*SUMTA**0. 3333 

IF  (EN1-0.0)  80,90,90 
C  IF  ANY  ACTIVATION  ENERGY  IS  NEGATIVE,  SET  IT  EQUAL  TO  ZERO 

80  EN1=0.0 
90  EN2=ENIN2 

EN2=ENIN2-0. 36* (0. 5+21 ./RDD)*EXP(-RDD/21 .)-3.6E-8*SUMTA**0. 3333 

IF  (EN2-0.0)  100,110,110 
100  EN2=0.0 
110  EN3=ENIN3 

EN3=ENIN3-0.36*(0.5+21./RDA)*EXP(-RDA/21 .)-3.6E-8*SUMTD**0.3333 

IF  (EN3-0.0)  120,130,130 
120  EN3=0.0 
130  EN4=ENIN4 

EN4=ENIN4-0.36*  (0.5+21 . /RDA)*EXP(-RDA/21 . )-3.6E-8*SUMTD**0. 3333 

IF  (EN4-0.0)  140,150,150 
140  EN4=0.0 
150  EN5=ENIN5 

EN5=ENIN5-0. 36*(0. 5+21 ./RDA)*EXP(-RDA/21 .)-3.6E-8*SUMTD**0. 3333 
IF  (EN5-0.0)  160,170,170 
160  EN5=0.0 

C  THE  FIRST  TEMPERATURE  WILL  BE  4.2  KELVIN 

170  T=4.2 

C  T  WILL  BE  INCREMENTED  IN  EQUAL  1000/T  INTERVALS 

RTEMP=1000./T 

C 

C  SET  UP  THE  OUTPUT  PAGE 

WRITE  (6,180) 

180       FORMAT  (1H1,32HELECTRICAL  PROPERTIES  OF  SILICON) 

IF  (NAME1 (1).EQ.6HC0NSTA)  GO  TO  320 
C  READ-OUT  ALL  PERTINENT  INPUT  DATA 

WRITE  (6,190) 
190       FORMAT  (1 H0,26HDOPANT  STATES  FOR  THIS  RUN) 

WRITE  (6,200) 

200       FORMAT  (1H0,77H        DOPANT  DENSITY  ENERGY (INPUT)  E 

2NERGY(CALC.)  DEGENERACY) 
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IF  (NAME1 (1) .EQ.6HDUMMY  )  60  TO  210 

WRITE  (6,250)  NAME1 ,DN1 ,ENIN1 ,EN1 ,DG1 
210      IF  (NAME2(1) .EQ.6HDUMMY  )  60  TO  220 

WRITE  (6,250)  NAME2,DN2,ENIN2,EN2,D62 
220      IF  (NAME3(1) .EQ.6HDUMMY  )  60  TO  230 

WRITE  (6,250)  NAME3,DN3,ENIN3,EN3,D63 
230      IF  (NAME4(1) .EQ.6HDUMMY  )  GO  TO  240 

WRITE  (6,250)  NAME4,DN4,ENIN4,EN4,DG4 
240      IF  (NAME5(1) .EQ.6HDUMMY  )  GO  TO  260 

WRITE  (6,250)  NAMES, DN5,ENIN5,EN5,D65 
250      FORMAT  (1X,2A6,6X,1PE9.3,7X,0PF7.4,8X,F8.4,10X,F8.4) 
260      WRITE  (6,270) 
270      FORMAT  (IHO) 

C  PROVIDE  COLUMN  HEADIN6S  FOR  OUTPUT  TABLE 

WRITE  (6,280) 

280      FORMAT  (1H0,129H      TEMP  1000/T  EF  N  P 

2  COND.  MOB  RH  RHO  D1   ION  D2  ION  A1 

3I0N        A2  ION) 
WRITE  (6,290) 
290      FORMAT  (1H  ) 

IF  (D0UBSW.NE.6H0N        )  GO  TO  350 
WRITE  (6,300) 

300      FORMAT  (1H  ,130H        EG  TEST  KT/Q  EG-EF  A3 

2I0N  SIGMA      MOB-(LAT)       +IONIZED       +NEUTRAL    MOB+(LAT)  +IONI 

3ZED  +NEUTRAL) 
WRITE  (6,310) 
310      FORMAT  (IHO) 
GO  TO  350 

C  COLUMN  HEADINGS  FOR  'CONSTANTS'  OPTION 

320      WRITE  (6,330) 

330      FORMAT  (1H0,128H      TEMP  1000/T  KT/Q  EG  DI 

2EL        M08-(LAT)     MOB+(LAT)       T**3/2  ME  MH  FN 

3C  FNV) 
WRITE  (6,340) 

340      FORMAT  (IHO) 

C  EFME  AND  EFMH  ARE  THE  EFFECTIVE  MASSES  OF  ELECTRONS  AND  HOLES 

C  1039  (1967).  EFME  AND  EFMH  EQS  FITTED  OVER  TEMP  RANGE  0  TO  600  K 

350      EFME= ( ( ( ( (4 . 54649E-1 7*T-9 . 66067E-1 4) *T+8 . 04032E-1 1 ) *T-3 . 3201 30E-8 ) 
2*T+6.83008E-6)*T-. 0001 61 708)*T+1. 0627 

352  EFMH=( ( (((1.11 997E-16*T-2. 596730E-13)*T+2.30049E-10)*T-9.67212E-8) 
2*T+1 .85678E-5)*T-.000523548)*T+. 590525 

353  CAYT=8.6173E-5*T 

C  CONST  =  2*(2*PI*K*(FREE  ELECTRON  MASS) /H**2)**(3/2) 

C0NST=4.829E15 
TP0WR=T**1.5 

C  FNC  AND  FNV  ARE  THE  DENSITY  OF  STATES  IN  THE  CONDUCTION  AND 

C  VALENCE  BANDS  RESPECTIVELY 

354  FNC=C0NST*(EFME**1 .5)*TP0WR 

355  FNV=C0NST*(EFMH**1 .5)*TP0WR 

C  EG  FROM  FIT  OF  DATA  FROM  MACFARLANE,  PHYS.   REV.  Ill,  1245  (1958) 

356  EG=(( (-3. 80977E-13*T+9.95E-10)*T-8. 701 10E-7)*T+. 0000323741 )*T+1. 15 
2556 

C  IN  'CONSTANTS'  OPTION,  WE  NEED  NOT  DETERMINE  EF 
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IF  (NAME1 (1) .EQ.6HC0NSTA)  60  TO  530 

EF=0.0 

EFG=1 .2 

RE=1 .OE-6 

AE=1 .OE-6 

IFLAG=0 

C  CALL  SUBROUTINE  TO  FIND  VALUE  OF  EF  WHICH  MAKES  TEST  VERY  SMALL 

CALL  ZEROIN  (TEST,EF,EFG,RE,AE,IFLAG) 
IF  (IFLAG.GT.2)  GO  TO  690 

C 

C  CALL  SUBROUTINE  TO  CALCULATE  MOBILITY 

CALL  MOB 

C  ERROR  AND/OR  WARNING  DIAGNOSTIC  MESSAGES 

IF  (ALPHA. LT. -2.)  GO  TO  370 

WRITE  (6,360)  ALPHA 
360      FORMAT  (20X,1 5H (EF-EG) /CAYT  =  ,F8.4, 

2  33H     ERROR  IN  N  MAY  EXCEED  1  PERCENT) 
370      IF  (BETA.LT.-2.)  GO  TO  390 

BETAN=-BETA 

WRITE   (6,380)  BETAN 
380      FORMAT  (20X,1 OHEF/CAYT  =  ,F8.4, 

2  33H     ERROR  IN  P  MAY  EXCEED  1  PERCENT) 
390      IF  (ISTAT.GE. 10000)  WRITE  (6,400) 

400       FORMAT  (20X,72HB0RN  APPROXIMATION  FOR  IONIZED  SCATTERING  NOT  VALID 
2  FOR  MAJORITY  CARRIER) 

C 

C  OUTPUT  THE  RESULTS  AND  SET  UP  THE  NEXT  RUN 

WRITE   (6,410)  T,RTEMP,EF 
410       FORMAT  ( F1 1 . 4,2X, F9 . 4, F1 1 . 6) 

IF  (CN.NE.O.)  WRITE  (6,420)  CN 

IF  (CP.NE.O.)  WRITE  (6,430)  CP 

IF   (UD.NE.O.)  WRITE  (6,440)  UD 

IF  (RH.NE.9.999999E19.AND.RH.NE.0.)  WRITE  (6,450)  RH 

IF  (RH.NE.9.999999E19.AND.RH0.NE.0.)  WRITE  (6,460)  RHO 

IF  (DI1.NE.0.)  WRITE  (6,470)  DI1 

IF  (0I2.NE.0.)  WRITE  (6,480)  DI2 

IF  (DI3.NE.0.)  WRITE  (6,490)  DI3 

IF  (DI4.NE.0.)  WRITE  (6,500)  DI4 
420  FORMAT  (1 H+,33X,1 PE1 1 . 3) 
430  FORMAT  (1H+,44X,1PE11 .3) 
440  FORMAT  (1 H+,55X,1 PEI 1 . 3) 
450  FORMAT  (1 H+,66X,1 PEI 1 .3) 
460  FORMAT  (1 H+,77X,1 PEI 1 . 3) 
470  FORMAT  (1H+,88X,1PE1 1 .3) 
480  FORMAT  (1 H+,99X,1 PEI 1 . 3) 
490  FORMAT  (1 H+,1 10X,1 PEI 1 .3) 
500      FORMAT  (1H+,121X,1PE10.3) 

IF  (D0UBSW.NE.6H0N        )  GO  TO  590 

0UT3=EG-EF 

TESTF=TEST(EF) 

WRITE  (6,510)  EG,TESTF,CAYT,0UT3 
510      FORMAT  (F11.6,1PE11.3,1PE11.4,2X,F9.6) 
IF  (DI5.NE.0.)  WRITE  (6,430)  DI5 
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IF  (RH.NE.9.999999E19.AND.SIGMA.NE.0.)  WRITE  (6,440)  SIGMA 
IF  (ULN.NE.O.)  WRITE  (6,450)  ULN 
IF  (0UT8.NE.0.)  WRITE  (6,460)  OUTS 
IF  (UDN.NE.O.)  WRITE   (6,470)  UON 
IF  (ULP.NE.O.)  WRITE  (6,480)  ULP 
IF  (0UT11 .NE.O.)  WRITE  (6,490)  0UT1 1 
IF  (UDP.NE.O.)  WRITE  (6,500)  UDP 
WRITE  (6,520) 
520      FORMAT  (1H0) 
GO  TO  590 

C  OUTPUT  INSTRUCTIONS  FOR  THE  'CONSTANTS'  OPTION 

530  RDIEL=11.4294*EXP(7.8E-5*T) 

IF  (T.GE.108.)  GO  TO  540 

ULN=1.89E7*T**(-1.52) 

GO  TO  550 
540  ULN=7.98E8*T**(-2.32) 
550      IF  (T.GE.72.)  GO  TO  560 

ULP=1.6E7*T**(-1.54) 

GO  TO  570 
560  ULP=2.3E9*T**(-2.7) 

570      WRITE  (6,580)  T,RTEMP,CAYT,EG,RDIEL,ULN,ULP,TPOWR,EFME,EFMH,FNC,FN 

2V 

580      FORMAT  (0P2F11 .4,1PE11 .4,2X,0PF9.4,F11 .4,1P7E11 .3) 
C  SET  UP  NEXT  TEMPERATURE 

C  THE  MAXIMUM  TEMPERATURE  OF  INTEREST  IS  500  DEGREES  KELVIN 

590      IF  (T-500.)  600,710,710 
600      IF  (RUNNO-2.)  610,620,630 
610  RUNN0=2. 

C  SECOND  TIME  AROUND,  THE  TEMPERATURE  WILL  BE  20  DEGREES  KELVIN 

T=20. 

RTEMP=1000./T 
GO  TO  350 

C  AFTER  THE  SECOND  TIME  AROUND,  1000. /T  IS  DECREMENTED  BY  2. 

620  RTEMP=RTEMP-2. 
T=1000./RTEMP 

C  AN  EVALUATION  AT  LIQUID  NITROGEN  TEMPERATURE  IS  DESIRED 

IF  (T.GE.77.)  GO  TO  640 

GO  TO  350 
630      IF  (T.EQ.77.)  GO  TO  650 

C  AN  EVALUATION  AT  ROOM  TEMPERATURE  (300  KELVIN)  IS  DESIRED 

IF  (RUNN0.EQ.4.)  GO  TO  680 

IF  (T.EQ.300.)  GO  TO  670 
C  FROM  77  KELVIN  TO  ROOM  TEMPERATURE,  1000/T  IS  DECREMENTED  BY  1. 

RTEMP=RTEMP-1. 

T=1000./RTEMP 

IF  (T.GT.300.)  GO  TO  660 

GO  TO  350 
640  T=77. 

RTEMP=1000./T 

RUNN0=3. 

GO  TO  350 
650  RTEMP=12. 

T=1000./RTEMP 
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GO  TO  350 
660  T=300. 

RTEMP=1000./T 

GO  TO  350 
670  RTEMP=3. 

T=1000./RTEMP 

RUNN0=4. 

GO  TO  350 

C  ABOVE  ROOM  TEMPERATURE     1000, /T  IS  DECREMENTED  BY  0.5 

680  RTEMP=RTEMP-0.5 

T=1000./RTEMP 

GO  TO  350 
690      WRITE  (6,700)  IFLAG 

700       FORMAT  (1  H0,51 HCHARGE  BALANCE  EQUATION  COULD  NOT  BE  SOLVED,  IFLAG= 
2,11) 
GO  TO  590 

710       IF   (D0UBSW.NE.6H0N        )  GO  TO  720 

C  INITIALIZE  CERTAIN  VARIABLES  FOR  NEXT  RUN 

ULN=0. 

0UT8=0. 

UDN=0. 

ULP=0. 

0UT11=0. 

UDP=0. 
720  CONTINUE 

IF   (NAME1 (1 ) .EQ.6HC0NSTA)  60  TO  10 

LOOPSW=0. 

HEADSW=6H0FF 

GO  TO  10 

C  ROUTINE  TO  SET-UP  OR  CANCEL  DOUBLE  OUTPUT  LINE  FORMAT 

C  DOUBSW  IS  A  SWITCH  USED  TO  CONTROL  THE  OUTPUT  FORMAT 

730  D0UBSW=6H0N 

IF  (NAME1  (2) .EQ.6H0N        )  GO  TO  10 

IF   (NAME1  (2) .EQ.6H  ON       )  GO  TO  10 
C  DEFAULT  VALUE  IS  TAKEN  TO  MEAN  'DOUBLE  OFF' 

D0UBSW=6H0FF 

GO  TO  10 

C  ROUTINE  TO  SET-UP  DOLOOP  OPTION 

740  DELDEN=DN1 
NDIT=ENIN1 

IF   (NAME1  (2) .EQ.6H  D1       )  L00PSW=1 . 

IF  (NAME1  (2).EQ.6H  D2       )  L00PSW=2. 

IF  (NAME1 (2) .EQ.6H  A1       )  L00PSW=3. 

IF   (NAME1  (2).EQ.6H  A2       )  L00PSW=4. 

IF  (NAME1 (2) .EQ.6H  A3       )  L00PSW=5. 

IF  (LOOPSW.EQ.O.)  WRITE  (6,750) 
750       FORMAT  (1H0,34HINVALID  ARGUMENT  OF  DOLOOP  COMMAND) 

GO  TO  10 
760  STOP 

END 
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FUNCTION  TEST  (EF) 

C 

C  THIS  EXTERNAL  FUNCTION  DEFINES  THE  PARAMETER  'TEST'. 

C  SUBROUTINE  ZEROIN  LOCATES  THE  VALUE  OF  EF  WHICH  MINIMIZES  'TEST'. 

DIMENSION  NAME1 (2) ,NAME2 (2) ,NAME3 (2) ,NAME4(2) , NAMES (2) 

COMMON  NAME1 ,NAME2,NAME3,NAME4, NAMES 

COMMON  DNI ,DN2,DN3,DN4,DNS,SUMTD,SUMTA 

COMMON  DI1,DI2,DI3,DI4,DIS 

COMMON  EN1,EN2,EN3,EN4,ENS 

COMMON  D61,DG2,DG3,DG4,DG5 

COMMON  T,CN, CP, EG, EFME,EFMH,FNC,FNV,CAYT,TPOWR,Q, ALPHA, BETA 
COMMON  UD,UDN,ULN,UDP,ULP, SIGMA, RHO,RH 
COMMON  0UT8,0UT11,ISTAT 

C 

C  COMPUTATION  OF  CN 

2  CN=0. 

4  ALPHA=(EF-EG)/CAYT 

IF  (ALPHA. GT.1.)  GO  TO  10 

IF  (ALPHA. LE. -80.)  60  TO  20 
C  THE  FOLLOWING  EQUATION  IS  USED  FOR  ALPHA  BETWEEN  -80  AND  +1 

8  CN=FNC/((EXP(-ALPHA))+0.27) 

GO  TO  20 

C  THE  FOLLOWING  EQUATION  IS  USED  FOR  ALPHA  GREATER  THAN  +1 

10        CN=FNC*0.7522S3*((ALPHA*ALPHA+1 .7)**0.7S) 

C  COMPUTATION  OF  CP 

20  CP=0. 

22  BETA=-EF/CAYT 

IF  (BETA. GT.1.)  GO  TO  30 

IF  (BETA. LE. -80.)  GO  TO  40 
C  THE  FOLLOWING  EQUATION  IS  USED  FOR  BETA  BETWEEN  -80  AND  +1 

26  CP=FNV/((EXP(-BETA))+0.27) 

60  TO  40 

C  THE  FOLLOWING  EQUATION  IS  USED  FOR  BETA  GREATER  THAN  +1 

30        CP=FNV*0.7S22S3*((BETA*BETA+1 .7)**0.7S) 
C  COMPUTATION  OF  DI1 

40  DI1=0. 

IF  (NAME1 (1).EQ.6HDUMMY  )  GO  TO  60 

P0WER=(EF-EG+EN1)/CAYT 

IF  (POWER. 6E. 80.)  60  TO  60 

IF  (POWER. LE. -80.)  60  TO  SO 

DI1=DN1/(1 .+(EXP(P0WER)*D61)) 

60  TO  60 
50  DI1=DN1 
C  COMPUTATION  OF  012 

60  DI2=0. 

IF  (NAME2(1) .EQ.6HDUMMY  )  60  TO  80 

P0WER=(EF-EG+EN2)/CAYT 

IF  (POWER. GE. 80.)  60  TO  80 

IF  (POWER. LE. -80. )  60  TO  70 

DI2=DN2/(1 .+(EXP(P0WER)*D62)) 

60  TO  80 
70  DI2=DN2 
C  COMPUTATION  OF  DI3 
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80  DI3=0. 

IF   (NAME3(1 ) .EQ.6HDUMMY  )  GO  TO  100 

P0WER=(EN3-EF)/CAYT 

IF  (POWER. GE. 80.)  GO  TO  100 

IF  (POWER. LE. -80.)  60  TO  90 

DI3=DN3/(1 .+(EXP(P0WER)*D63)) 

GO  TO  100 
90  DI3=DN3 
C  COMPUTATION  OF  DI4 

100  DI4=0. 

IF   (NAME4(1 ) .EQ.6HDUMMY  )  GO  TO  120 

P0WER=(EN4-EF)/CAYT 

IF  (POWER. GT. 80.)  GO  TO  120 

IF  (POWER. LE. -80. )  GO  TO  110 

DI4=DN4/(1 .+(EXP(P0WER)*DG4)) 

GO  TO  120 
110  DI4=DN4 
C  COMPUTATION  OF  DI5 

120  DI5=0. 

IF   (NAME5(1 ) .EQ.6HDUMMY  )  GO  TO  140 

P0WER=(EN5-EF)/CAYT 

IF   (POWER. GE. 80.)  GO  TO  140 

IF  (POWER. LE. -80.)  GO  TO  130 

DI5=DN5/(1  .  +  (EXP(P0WER)*D65)) 

GO  TO  140 
130  DI5=DN5 

C  TESTING  TO  SEE  HOW  CLOSE  TO  NEUTRALITY  WE  HAVE  COME 

C  THE  PARAMETER  TEST  IS  THE  NET  DENSITY  OF  RESIDUAL  CHARGES 

C  WHEN  THE  FERMI  LEVEL  IS  AT  ITS  PRESENT  POSITION 

140  TEST=(CN+DI3+DI4+DI5)-(CP+DI1+DI2) 

RETURN 

END 
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SUBROUTINE  ZEROIN(F,B,C,RE,AE,IFLAG) 

C 

C  SANDIA  MATHEMATICAL  PROGRAM  LIBRARY 

C  MATHEMATICAL  COMPUTING  SERVICES  DIVISION  5422 

C  SANDIA  LABORATORIES 

C  P.  0.  BOX  5800 

C  ALBUQUERQUE,  NEW  MEXICO  87115 

C  CONTROL  DATA  6600  VERSION  4.5,  1  NOVEMBER  1971 

C 

C  MODIFIED  TO  RUN  AT  NBS  BY  D.  KAHANER  DIV  713 
C 

C  ABSTRACT 

C  ZEROIN  SEARCHES  FOR  A  ZERO  OF  A  FUNCTION  F(X)  BETWEEN 

C  THE  GIVEN  VALUES  B  AND  C  UNTIL  THE  WIDTH  OF  THE  INTERVAL 

C  (B,C)  HAS  COLLAPSED  TO  WITHIN  A  TOLERANCE  SPECIFIED  BY 

C  THE  STOPPING  CRITERION,  ABS(B-C)   .LE.  2.*(RW*ABS(B)+AE) . 
C 

C  DESCRIPTION  OF  ARGUMENTS 

C  F          -  NAME  OF  THE  REAL  VALUED  EXTERNAL  FUNCTION.     THIS  NAME 

C  MUST  BE  IN  AN  EXTERNAL  STATEMENT  IN  THE  CALLING 

C  PROGRAM.     F  MUST  BE  A  FUNCTION  OF  ONE  REAL  ARGUMENT. 

C  B          -  ONE  END  OF  THE  INTERVAL  (B,C) .     THE  VALUE  RETURNED  FOR 

C  B  USUALLY  IS  THE  BETTER  APPROXIMATION  TO  A  ZERO  OF  F. 

C  C          -  THE  OTHER  END  OF  THE  INTERVAL  (B,C) 

C  RE        -  RELATIVE  ERROR  USED  FOR  RW  IN  THE  STOPPING  CRITERION. 

C  IF  THE  REQUESTED  RE  IS  LESS  THAN  MACHINE  PRECISION, 

C  THEN  RW  IS  SET  TO  APPROXIMATELY  MACHINE  PRECISION. 

C  AE        -  ABSOLUTE  ERROR  USED  IN  THE  STOPPING  CRITERION.     IF  THE 

C  GIVEN  INTERVAL  (B,C)  CONTAINS  THE  ORIGIN,  THEN  A 

C  NONZERO  VALUE  SHOULD  BE  CHOSEN  FOR  AE. 

C  IFLAG  -  RETURNS  A  STATUS  OF  THE  RESULTS  INDICATING  WHICH 

C  OF  THE  FOLLOWING  CONDITIONS  HOLD. 

C  A  -  ABS(B-C)   .LE.  2.*(RW*ABS CB)+AE) 

C  B  -  F(B)  *  FCC)   .LT.  0. 

C  C  -  ABSCFCB))   .LE.  ABSCFCO) 

C  D  -  ABSCFCB      ))   .LE.  MAXCABSCFCB     )),ABSCFCC  ))) 

C  OUT  IN  IN 

C  E  -  NUMBER  OF  EVALUATIONS  OF  FCX)   .LE.  500 

C  =1   INDICATES  NORMAL  CASE.     ALL  CONDITIONS  ABOVE  HOLD. 

C  =2  INDICATES  FCB)  =  0.     CONDITION  A  MAY  NOT  HOLD. 

C  =3  INDICATES  CONDITIONS  A,  B,  C,  AND  E  HOLD  BUT  D  DOES 

C  NOT.     CB,C)  PROBABLY  CONTAINS  A  SINGULAR  POINT  OF  F. 

C  =4  INDICATES  CONDITIONS  A  AND  E  HOLD  BUT  B  DOES  NOT. 

C  A  LOCAL  MINIMUM  OF  FCX)  I.M  CB,C)  MAY  HAVE  BEEN  FOUND 

C  =5  INDICATES  SEARCH  WAS  ABORTED  WHEN  CONDITION  E  FAILED 

C 

C  REFERENCES 

C  1.     L  F  SHAMPINE  AND  H  A  WATTS,  ZEROIN,  A  ROOT-SOLVING  CODE, 

C  SC-TM-70-631,  SEPT  1970. 

C  2.     T  J  DEKKER,  FINDING  A  ZERO  BY  MEANS  OF  SUCCESSIVE  LINEAR 

C  INTERPOLATION,  *CONSTRUCTIVE  ASPECTS  OF  THE  FUNDAMENTAL 

C  THEOREM  OF  ALGEBRA*,  EDITED  BY  B  DEJON  AND  P  HENRICI,  1969. 

C  3.     L  F  SHAMPINE  AND  R  C  ALLEN,  NUMERICAL  COMPUTING— 
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AN  INTRODUCTION,  1973. 


INITIALIZE  FOR  UNIVAC  1108 
DATA  ER/6.0E-8/ 
INITIALIZE  IN  PORTABLE  MANNER 
DATA  ER/0.0/ 

IF  (ER  .EQ.  0.0)  ER  =  4.  *  (R1MACH  (4)  ) 

RW=AMAX1 (RE,ER) 

IC=0 

ACBS=ABS(B-C) 
A=C 

FA=F(A) 
FB=F(B) 
FC  =  FA 
K0UNT=2 

FX=AMAX1 (ABS(FB),ABS(FC)) 

IF  (ABS(FC)   .GE.  ABS(FB))  GO  TO  2 

PERFORM  INTERCHANGE 

A=B 

FA=FB 

B=C 

FB=FC 

C=A 

FC  =  FA 

CMB=0.5*(C-B) 
ACMB=ABS(CMB) 
TOL=RW*ABS(B)+AE 

TEST  STOPPING  CRITERION 

IF  (ACMB  .LE.  TOL)  GO  TO  10 

CALCULATE  NEW  ITERATE  IMPLICITLY  AS  B+P/Q 
WHERE  WE  ARRANGE  P  .GE.  0. 

THE  IMPLICIT  FORM  IS  USED  TO  PREVENT  OVERFLOW. 

P=(B-A)*FB 

Q=FA-FB 

IF  (P  .GE.  0.)  GO  TO  3 

P=-P 

Q=-Q 

UPDATE  A  AND  CHECK  FOR  SATISFACTORY  REDUCTION 

IN  THE  SIZE  OF  OUR  BOUNDING  INTERVAL. 

A=B 

FA=FB 

IC=IC+1 

IF  (IC   .LT.  4)  GO  TO  4 

IF   (8.*ACMB  .GE.  ACBS)  GO  TO  6 

IC=0 

ACBS=ACMB 
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C  TEST  FOR  TOO  SMALL  A  CHANGE 

4  IF  (P  .GT.  ABS(Q)*TOL)  60  TO  5 

C 

C  INCREMENT  BY  TOLERANCE 

B=B+$IGN(TOL,CMB) 
60  TO  7 

C 

C  ROOT  0U6HT  TO  BE  BETWEEN  B  AND  (C+B)/2. 

5  IF  (P  .6E.  CMB*Q)  60  TO  6 

C 

C  INTERPOLATE 
B=B+P/Q 
60  TO  7 

C 

6  B=0.5*(C+8) 
C  BISECT 

C 

C  HAVE  COMPLETED  COMPUTATION  FOR  NEW  ITERATE  B 

7  FB=F(B) 

IF  (FB  .EQ.  0.)  60  TO  11 
K0UNT=K0UNT+1 

IF  (KOUNT  .6T.  500)  60  TO  15 

C 

C  DECIDE  WHETHER  NEXT  STEP  IS  INTERPOLATION  OR  EXTRAPOLATION 

IF  (SI6N(1 .0,FB)   .NE.  SI6N (1 .0, FC) )  60  TO  1 
C=A 
FC  =  FA 
60  TO  1 

C 
C 

C  FINISHED.  PROCESS  RESULTS  FOR  PROPER  SETTIN6  OF  IFLA6 

C 

10  IF  (FB*FC  .6T.  0.)  60  TO  13 
IF  (ABS(FB)   .6T.   FX)  60  TO  12 
IFLA6  =  1 

RETURN 

11  IFLA6  =  2 
RETURN 

12  IFLA6  =  3 
RETURN 

13  IFLA6  =  4  I 
RETURN 

15  IFLA6  =  5 
RETURN 
END 
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SUBROUTINE  MOB 

C 

C  SUBROUTINE  TO  EVALUATE  CARRIER  MOBILITIES  AND  RELATED  ITEMS 

C  THIS  SUBROUTINE  CALLS  SUBROUTINE  SICI  (SI,CI,X) 

DIMENSION  NAME1 (2) ,NAME2 (2) ,NAME3(2) ,NAME4(2) ,NAME5 (2) 

COMMON  NAME1,NAME2,NAME3,NAME4, NAMES 

COMMON  DN1,DN2,DN3,DN4,DN5,SUMTD,SUMTA 

COMMON  DI1,DI2,DI3,DI4,DI5 

COMMON  EN1,EN2,EN3,EN4,EN5 

COMMON  DG1,DG2,DG3,DG4,DG5 

COMMON  T,CN, CP, EG, EFME,EFMH,FNC,FNV,CAYT,TPOWR,Q, ALPHA, BETA 

COMMON  UD,UDN,ULN,UDP,ULP,SIGMA,RHO,RH 

COMMON  0UT8,0UT11,ISTAT 
C  UPON  RETURNING,  ISTAT  INDICATES  THE  STATUS  OF  THE  COMPUTATIONS 

C  EACH  DIGIT  (POSITION)  OF  ISTAT  IS  USED  INDIVIDUALLY 

C 

ISTAT=0 

C  DIEL  =  4*PI*(ABS0LUTE  DIELECTRIC  CONSTANT  IN  FARADS/METER) 

C  RELATIVE  DIELECTRIC  CONSTANT  AT  300  KELVIN  =  11.7 

C  DUNLAP  AND  WATTERS,  PHYS.REV.  92,  1396-1397  (1953) 

C  TEMPERATURE  DEPENDENCE  CALCULATED  FROM  INDEX  OF  REFRACTION  DATA 

C  CARDONA  ET  AL,  J.   PHYS.   CHEM.   SOLIDS.  8,  204-206  (1959) 

2  DIEL=1 .2711E-9*EXP(7.8E-5*T) 

C 

C  SEE  PAPERS  BY  LI  AND  THURBER,  SSE  20,  609  AND  LI,  SSE  21,  1109 

SUMID=DI1+DI2 

SUMIA=DI3+DI4+0I5 

CI0N=DI1+DI2+DI3+DI4+DI5 

CNUT=DN1+DN2+DN3+DN4+DN5-CI0N 
C  WE  MUST  GUARD  AGAINST  ZERO  DIVISORS  IN  THE  FOLLOWING  COMPUTATIONS 

C  SOME  OF  THESE  ZERO  DEVISORS  WILL  BE  FLAGGED  BY  SETTING  ISTAT 

IF  (CION.LE.O.)  CI0N=1. 

IF   (CNUT.LE.O.)  CNUT=1. 
C  CALCULATION  OF  ELECTRON  MOBILITY 

C  LATTICE  SCATTERING  (WITH  CORRECTION  FOR  ELECT . -ELECT.  SCATTERING) 

C  EQ  FOR  CORRECTION  DEPENDS  ON  DONOR  DENSITY 

C  EQS  FOR  TEMPERATURE  DEPENDENCE  OF  ULN  FROM  NORTON,  ET  AL. 

C  SEE  PAPER  IN  PHYS.  REV.  B8,  5632  (1973) 

C  VALUE  OF  ULN  AT  300  K  FROM  THURBER,  ET  AL.  TO  BE  PUBLISHED  IN  J  ECS 

IF  (SUMTD.GE.3.E17)  60  TO  20 

IF  (T.GE.108.)  GO  TO  10 
8  ULN=1 .89E7*T**(-1 .52)*(1 .-0.08*SUMTD/2.E17) 

GO  TO  40 

10        ULN=7.98E8*T**(-2.32)*(1 .-0.08*SUMTD/2.E17) 

GO  TO  40 
20        IF  (T.6E.108.)  60  TO  30 
22        ULN=1 .89E7*(T**(-1 .52))*0.884 

GO  TO  40 

30  ULN=7.98E8*(T**(-2.32))*0.884 

C  ADD  ON  THE  EFFECTS  OF  IONIZED  IMPURITY  SCATTERIN6 

C  WE  MUST  ALLOW  FOR  THE  FACT  THAT  SUMTD  MIGHT  VANISH  IN  CNSTAR  EQ. 

40        IF  (SUMTD. LE.O.)  CNSTAR=CN+CP 
IF  (SUMTD. LE.O.)  60  TO  50 
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42        CNSTAR=CN+CP+(1 .-SUMID/SUMTD)*SUMID 

C  BN  IS  LOG  TERM  IN  MODIFIED  BROOKS-HERRING  IONIZED  IMPURITY  MOB.  EQ 

C  BN  =  6*DIEL*MASS*(1 .3806E-23)**2*T*T/(PI*Q*Q*HBAR*HBAR*CNSTAR) 

C  WE  MUST  ALLOW  FOR  THE  FACT  THAT  CNSTAR  MIGHT  VANISH  IN  EQ  FOR  BN 

C  BN  IS  SET  TO  A  LARGE  NUMBER  AND  FIRST  POSITION  OF  ISTAT 

C  IS  SET  TO  UNITY 

50        IF  (CNSTAR. LE.O.)  BN=9.999999E1 9 

IF  (CNSTAR. LE.O.)  ISTAT=ISTAT+1 

IF  (CNSTAR. LE.O.)  GO  TO  80 
52        BN=1 .14E24*DIEL*T*T/CNSTAR 

C  IF  BN  IS  SMALL,  BORN  APPROXIMATION  IS  NOT  VALID.     IN  THIS  CASE, 

C  FIFTH  POSITION  OF  ISTAT  IS  INCREMENTED. 

IF  (BN-10.)  60,60,80 
60  ISTAT=ISTAT+10000 

IF  (CION-1.0)  70,80,80 
70  UDN=ULN 

0UT8=UDN 

GO  TO  110 
80        GN=AL0G(8N+1 .)-BN/(BN+1 .) 

C  CALCULATE  GAMMA  FOR  EFFECT  OF  E-E  SCATTERING  ON  ION.   IMP.  MOBILITY 

IF  (CNSTAR. LE.O.)  GAMMA=0.632 

IF  (CNSTAR. LE.O.)  GO  TO  90 

IF  (SUMID.LE.O.)  GAMMA=0.432 
86        IF  (SUMID.GT.O.)  GAMMA=(SUMID/CNSTAR)*(1 .-EXP(-CNSTAR/SUMID) ) 

IF  (GAMMA. LT. 0.432)  GAMMA=0.432 
90  UIN=GAMMA*7.3E17*TP0WR/(CI0N*GN) 

IF  (UIN.LE.O.)  RATI0N=1.E6 

IF  (UIN.LE.O.)  60  TO  100 
98  RATI0N=6.*ULN/UIN 
100      XN=SQRT (RATION) 

C  COMBINE  LATTICE  AND  IONIZED  IMPURITY  SCATTERING  BY  THE  MIXED- 

C  SCATTERING  FORMULA  OF  DEBYE  AND  CONWELL,  PHYS.   REV.  93, 

C  698-706  (1954).  SUBROUTINE  SICI  PERFORMS  THIS  FUNCTION 

CALL  SICI  (SI,CI,XN) 
104      UDN=ULN*(1 .0+RATI0N*(CI*C0S(XN)+SIN(XN)*(SI-1 .5707963))) 

0UT8=UDN 

C  ADD  ON  THE  EFFECTS  OF  NEUTRAL  IMPURITY  SCATTERING 

C  NEUTRAL  IMPURITY  MOBILITY  AS  PER  ER6INS0Y,  PHYS  REV  79,1013  (1950) 

C  UNIS  =  0.02*PI**3*Q**3*MASS/(5*DIEL*CNUT*H**3) 

C  UNIS  CONTAINS  GEOMETRIC  MEAN  MASS  WITH  VALUE  OF  £FME/6**0. 667 

C  SEE  BROOKS,  ADV.   IN  ELECTRONICS  AND  ELECTRON  PHYSICS  7,  161  (1955) 

110  UNIS=4.85E11*EFME/(DIEL*CNUT) 

C  SCLAR'S  CORRECTION  TO  UNIS,  SEE  PHYS  REV  104,  1559  (1956) 

C  SEE  ALSO  PAPER  BY  LI  AND  THURBER,  SSE  20,  609  (1977) 

112      EI=0.215*EFME/((DIEL/1 .1121E-10)**2*CAYT) 
114      UNIS=UNIS*0.82*(2*SQRT(1/EI) /3+SQRT(EI) /3. ) 
116       UDN=1/(1 /UDN+1/UNIS) 

C  THAT  COMPLETES  THE  COMPUTATION  OF  THE  ELECTRON  MOBILITY 

C  CALCULATION  OF  HOLE  MOBILITY 

C  LATTICE  SCATTERING  WITH  CORRECTION  FOR  HOLE-HOLE  SCATTERING 

C  EQ  FOR  ULP  FOR  T.LT.72  FROM  BRAGGINS,  PHD  THESIS, 

C  SYRACUSE  UNIVERSITY,  1975 

C  EQ  FOR  ULP  FOR  T.GE.72  FROM  LUDWIG  AND  WATTERS, 


32 


C  PHYS  REV  101,  1699  (1956) 

IF  (SUMTA.GE.3.E17)  GO  TO  130 

IF  (T.GE.72.)  GO  TO  120 
118      ULP=1.6E7*T**(-1 .54)*(1 .-0.08*SUMTA/2.E17) 

GO  TO  150 

120      ULP=2.3E9*T**(-2.7)*(1 .-0.08*SUMTA/2. E17) 

GO  TO  150 
130      IF  (T.GE.72.)  GO  TO  140 
132  ULP=1.6E7*(T**(-1.54))*0.884 

GO  TO  150 
140  ULP=2.3E9*(T**(-2.7))*0.884 


C  ADD  ON  THE  EFFECTS  OF  IONIZED  IMPURITY  SCATTERING 

C  WE  MUST  ALLOW  FOR  THE  FACT  THAT  SUMTA  MIGHT  VANISH  IN  CPSTAR  EQ. 

150  IF  (SUMTA. LE.O.)  CPSTAR=CN+CP 

IF  (SUMTA. LE.O.)  GO  TO  160 

152  CPSTAR=CP+CN+(1 .-SUMIA/SUMTA)*SUMIA 

C  BP  IS  LOG  TERM  IN  BROOKS-HERRING  IONIZED  IMPURITY  MOBILITY 

C  BP  =  6*DIEL*MASS*(1 .3806E-23)**2*T*T/(PI*Q*Q*HBAR*HBAR*CPSTAR) 

C  MASS  FOR  BP  OF  EACH  BAND  FROM  LI,  NBS  SPECIAL  PUBLICATION  400-47 

C  NO  TEMPERATURE  DEPENDENCE  FOR  MASS,  300  K  VALUE  USED 

C  IS  SET  TO  1 

C  WE  MUST  ALLOW  FOR  THE  FACT  THAT  CPSTAR  MIGHT  VANISH  IN  EQ  FOR  BP 
C  BP  IS  SET  TO  A  LARGE  NUMBER  AND  THIRD  POSITION  OF  ISTAT 


160      IF  (CPSTAR. LE.O.)  BP1 =9 . 999999E1 9 

IF  (CPSTAR. LE.O.)  BP2=9. 999999E1 9 

IF  (CPSTAR. LE.O.)  BP3=9 . 999999E1 9 

IF  (CPSTAR. LE.O.)   ISTAT=ISTAT+1 00 

IF  (CPSTAR. LE.O.)  GO  TO  190 
162      BP1=1 .16E24*0.568*DIEL*T*T/CPSTAR 

BP2=1 . 1 6E24*0 . 41 2*DIEL*T*T/ CPSTAR 

BP3=1 .16E24*0.246*DIEL*T*T/CPSTAR 
C  IF  BP  IS  SMALL,  BORN  APPROXIMATION  IS  NOT  VALID.     IN  THIS  CASE, 

C  FIFTH  POSITION  OF  ISTAT  IS  INCREMENTED. 

IF  (BP1-10.)  170,170,190 
170  ISTAT=ISTAT+10000 

IF  (CION-1.0)  180,190,190 
180  UDP=ULP 

0UT11=UDP 

GO  TO  220 
1  90      GP1  =ALOG  (BP1  +1 . )  -BP1  /  (BP1  +1 .  ) 

GP2=AL0G (BP2+1 . ) -BP2/ (BP2+1 . ) 

GP3=AL0G (BP3+1 . )-BP3/ (BP3+1 .  ) 
C  UIP  =  0.01*(2**3.5)*DIEL*DIEL*SQRT(MD)*(1 .38E-23*T)**1 .5/ 

C  ((PI**1 .5)*Q*Q*Q*MC*CI0N*GP) 

192      UIP1=2.66E35*DIEL*DIEL*TP0WR*SQRT(0.568) /(0.448*CI0N*GP1) 
UIP2=2.66E35*DIEL*DIEL*TP0WR*SQRT(0.412)/(0.532*CI0N*GP2) 
UIP3=2.66E35*DIEL*DIEL*TP0WR*SQRT (0.246) /(0.253*CI0N*GP3) 
C  WEIGHT  INDIVIDUAL  MOBILITIES  BY  HOLE  DENSITY  IN  THAT  BAND 

196  UIP= (UIP1+UIP2* (0.41 2/0. 568) **1.5+UIP3*EXP(-0.044/CAYT)* (0.246/0. 5 
268)**1 .5)/(1 . +(0. 412/0. 568)**1 .5+EXP(-0. 044/CAYT)*(0. 246/0. 568)**1 
3.5) 

C  CALCULATE  GAMMA  FOR  EFFECT  OF  H-H  SCATTERING  ON  ION.   IMP.  MOBILITY 

IF  (CPSTAR. LE.O.)  GAMMA=0.632 
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IF  (CPSTAR.LE.O.)  60  TO  200 

IF  (SUMIA.LE.O.)  GAMMA=0.432 
198       IF  (SUMIA.GT.O.)  6AMMA=(SUMIA/CPSTAR)*(1 .-EXP(-CPSTAR/SUMIA)) 

IF  (GAMMA. LT. 0.432)  GAMMA=0.432 
200  UIP=GAMMA*UIP 

IF  (UIP.LE.O.)  RATI0P=1.E6 

IF  (UIP.LE.O.)  GO  TO  210 
208  RATI0P=6.*ULP/UIP 
210  XP=SQRT(RATIOP) 


C  COMBINE  LATTICE  AND  IONIZED  IMPURITY  SCATTERING  BY  THE  MIXED- 

C  SCATTERING  FORMULA  OF  DEBYE  AND  CONWELL,  PHYS.   REV.  93, 

C  698-706  (1954).  SUBROUTINE  SICI  PERFORMS  THIS  FUNCTION 

CALL  SICI  (SI,CI,XP) 
214      UDP=ULP*(1 .0+RATI0P*(CI*C0S(XP)+SIN(XP)*(SI-1 .5707963))) 

0UT11=UDP 

C  ADD  ON  THE  EFFECTS  OF  NEUTRAL  IMPURITY  SCATTERING 

C  NEUTRAL  IMPURITY  MOBILITY  AS  PER  ERGINSOY,  PHYS  REV  79,1013  (1950) 

C  UNIS  =  0.02*PI**3*Q**3*MASS/(5*DIEL*CNUT*H**3) 

C  VALENCE  BAND  OF  SILICON  HAS  THREE  SEPARATE  COMPONENTS 

C  PROCEDURE  TAKEN  FROM  LI,  NBS  SPECIAL  PUBLICATION  400-47 

C  FOR  THIS  CALCULATION,  MASS  VALUES  OF  LI  AT  300  K  ARE  USED 

220  UNIS1=1 .60E12*0.568/(DIEL*CNUT) 
UNIS2=1.60E12*0.412/(DIEL*CNUT) 
UNIS3=1 .60E12*0.246/(DIEL*CNUT) 

C  SCLAR'S  CORRECTION  TO  UNIS,  SEE  PHYS  REV  104,  1559,  (1956) 

221  EI1=0.71*0.568/((DIEL/1 .1121E-10)**2*CAYT) 
EI2=0.71*0.412/((DIEL/1 .1121E-10)**2*CAYT) 
EI3=0.71 *0.246/((DIEL/1.112lE-10)**2*CAYT) 

222  UNIS1=UNIS1*0.82*(2.*SQRT(1 ./EI1)/3.+SQRT(EI1) /3.) 
UNIS2=UNIS2*0.82*(2.*SQRT(1 ./EI2)/3.+SQRT(EI2)/3.) 
UNIS3=UNIS3*0.82*(2.*SQRT(1 ./EI3) /3. +SQRT (EI3) /3. ) 

C  WEIGHT  INDIVIDUAL  MOBILITIES  BY  HOLE  DENSITY  IN  THAT  BAND 

224      UNIS=(UNIS1+UNIS2*  (0.41 2 /0.568)**1.5+UNIS3*EXP(-0.044/CAYT)* (0.246 


2/0. 568)**1 .5)/(1. +(0.412/0. 568)**1 .5+EXP(-0.044/CAYT)*(0. 246/0. 568 
3)*1.5) 
226      UDP=1 /(I /UDP+1 /UNIS) 

C  THAT  COMPLETES  THE  COMPUTATION  OF  THE  HOLE  MOBILITY 

IF  (CN.LE.O. .AND.CP.LE.O. )  GO  TO  230 

RH=((CP*UDP*UDP-CN*UDN*UDN)/(Q*(CP*UDP+CN*UDN)))/(CP*UDP+CN*UDN) 
SIGMA =Q*(CN*UDN+CP*UDP) 
RH0=1 ./SIGMA 

UD=(CN*UDN+CP*UDP)/(CN+CP) 
C  END  OF  THE  MOBILITY  AND  RELATED  COMPUTATIONS 

C 

C  UPON  RETURNING,   ISTAT  INDICATES  THE  STATUS  OF  THE  COMPUTATIONS 

C  FIRST  DIGIT  UNITY,  CNSTAR  WAS  EQUAL  TO  ZERO 

C  THIRD  DIGIT  UNITY,  CPSTAR  WAS  EQUAL  TO  ZERO 

C  FIFTH  DIGIT  1  OR  2,  BORN  APPROXIMATION  NOT  VALID 

C  9.999999E19  IS  A     LARGE  NUMBER  INTRODUCED  TO  AVOID  ZERO  DIVISORS 

C 

230  RETURN 
END 
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SUBROUTINE  SICI (SI,CI,X) 

C 

C  THIS  ROUTINE  COMBINES  LATTICE  AND  IONIZED  MOBILITY  SCATTERING 

C  ACCORDING  TO  THE  MIXED-SCATTERING  FORMULA  OF  DEBYE  AND 

C  CONWELL,  PHYS.  REV.  93,  693-706  (1954) 

C 

C  TEST  ARGUMENT  RANGE 

Z=ABS(X) 

IF  (Z-4.)  10,10,50 
C  Z  IS  NOT  GREATER  THAN  4. 

10  Y=Z*Z 

S I=X* (((((( . 979421 54E-1 1 *Y-. 22232633E-8)*Y+. 30561 233E-6)*Y-. 283414 
260E-4)*Y+.16666582E-2)*Y-.55555547E-1)*Y+1.) 
C  TEST  FOR  LOGARITHMIC  SINGULARITY 

IF  (Z)  30,20,30 
20  CI=-1.E38 

RETURN 

30        CI=0. 57721 566+ALOG (Z)-Y*( ( ( ( (-.13869851 E-9*Y+.26945842E-7)*Y-.3095 

22207E-5 ) * Y+ . 231 46303E-3) *Y- . 1 041 6642E-1 ) *Y+ . 24999999) 
40  RETURN 

C  Z  IS  GREATER  THAN  4. 

50  SI=SIN(Z) 

Y=COS(Z) 

Z=4./Z 

U=((((((((.40480690E-2*Z-.022791426)*Z+.055150700)*Z-. 07261 641 8)*Z 
2+. 0498771 59)*Z-. 333251 86E-2)*Z-. 0231461 68)*Z-.11349579E-4)*Z+. 0625 
300111)*Z+.25839886E-9 

V=( (((((( ( (-. 0051 086993*Z+. 0281 91 786)*Z-. 065372834) *Z+. 079020335)* 
2Z-. 0440041 55 )*Z-.0079455563)*Z+. 02601 2930)*Z-.37640003E-3)*Z-. 0312 
3241 78) *Z-.66464406E-6)*Z+. 25000000 
CI=Z*(SI*V-Y*U) 
SI=-Z*(SI*U+Y*V)+1 .5707963 
C  TEST  FOR  NEGATIVE  ARGUMENT 

IF  (X)  60,40,40 
C  X  IS  LESS  THAN  -4. 

60  SI=-SI 
RETURN 
END 
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5.     Illustrative  Output  Data  Cards  and  Examples 
5.1     Data  Cards  Used  to  Generate  Output  Listings 

The  output  examples  to  follow  illustrate  a  no-option  run,   the  DOUBLE  ON 
format,   the  DOLOOP  XX  option,  and  the  CONSTANTS  option.     Following  is  a 
listing  of  the  complete  deck  of  dopant  state  data  and  control  cards  that 
used  to  generate  these  output  examples.     Lines  separate  the  card  sets  for 
each  example. 


PHOSPHORUS 

1 

.00000E15 

0. 

045 

2. 

DUMMY 

DUMMY 

BORON 

1 

.00000E13 

0. 

045 

4. 

DUMMY 

DOUBLE  ON 

PHOSPHORUS 

1 

.00000E18 

0. 

045 

2. 

DUMMY 

DUMMY 

BORON 

1 

.00000E13 

0. 

045 

4. 

DUMMY 

DOUBLE  OFF 

DOLOOP  A2 

5 

.OOOOOEU 

4. 

PHOSPHORUS 

1 

.00000E15 

0. 

045 

2. 

DUMMY 
DUMMY 

BORON  0.045  4. 

DUMMY  

CONSTANTS  

STOP  
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6.  CONCLUSIONS 


A  FORTRAN  program  has  been  presented  for  calculating  the  electrical  proper- 
ties of  silicon  containing  five  dopant  states.     The  charge-balance  equation 
is  solved  to  arrive  at  the  Fermi  energy  which  determines  the  carrier  densi- 
ties.    An  effective  conductivity  mobility  is  calculated  for  both  electrons 
and  holes,   taking  into  account  the  significant  scattering  mechanisms,  and 
then  the  resistivity  is  computed  from  the  mobilities  and  carrier  densities. 
An  attempt  has  been  made  to  include  in  the  algorithm  the  temperature  and  dop- 
ant density  dependence  of  as  many  parameters  as  possible.     In  some  cases 
(e.g.,   the  temperature  dependence  of  the  indirect  energy  gap  of  silicon), 
there  are  numerous  experimental  results  in  the  literature   [27] .     In  cases 
where  there  is  significant  disagreement  in  the  literature  (e.g,,   the  varia- 
tion of  dopant  activation  energy  with  dopant  density),   there  is  room  for  im- 
provement.    In  a  few  cases  where  little  or  nothing  is  known  (e.g.,   the  varia- 
tion of  the  Hall  scattering  factor  with  temperature  and  dopant  density),  no 
attempt  has  been  made  to  include  these  effects.     Therefore,   the  user  of  this 
program  is  urged  to  review  the  literature  regularly  for  additional  or  more 
comprehensive  measurements  that  would  improve  the  accuracy  of  the  computed 
results.     The  authors  invite  any  suggestions  and  comments  in  this  regard. 
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